1、测序质量评估
## conda安装fastqc
fastqc *.fastq.gz -t 8 -o fastqc_out/
可选(合并质控报告)
## conda安装multiqc (可能会要求安装latex,为生成PDF文件所需)
multiqc .
2、数据过滤(重复、接头等)
## conda安装fastp(fastp也会生成过滤前后的质量报告,如果需要的话,需要对报告文件重命名)
fastp -i in.R1.fq.gz -I in.R2.fq.gz -o out.R1.fq.gz -O out.R2.fq.gz
3、下载参考基因组(以玉米v5参考基因组为例)
## 下载参考基因组数据
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/902/167/145/GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0/GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna.gz
##解压基因组文件
gunzip GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna.gz
##更改名字
mv GCF_902167145.1_Zm-B73-REFERENCE-NAM-5.0_genomic.fna genome.fasta
4、构建索引
## 构建BWA索引
bwa index genome.fasta
bwa index的用法和参数:
bwa index [-p prefix] [-a algoType] <in.db.fasta>
-p #输出数据库的前缀,默认与输入文件名一致
-a #设置参考基因组构建index的算法(is和bwtsw两种)
bwa构建参考基因组index的两种算法:
-is 用于构建后缀数组的线性时间算法:缺点:比较耗内存(需要参考基因组5.37倍的内存);优点:速度快。
is是构建index的默认算法,但是不适合 “参考基因组” 大于2GB的文件。
-bwtsw 使用BWT-SW算法,多用于人类基因组。
##构建基因组索引,该命令会在运行目录下创建一个genome.fai索引文件
samtools faidx genome.fasta
##构建字典序文件,生成dict文件
gatk CreateSequenceDictionary -R example.fasta -O example.dict 或
picard CreateSequenceDictionary -R genome.fasta
GATK (全称The Genome Analysis Toolkit)是Broad Institute开发的二代重测序数据分析软件,是基因分析的工具集。在4.0以后,GATK包含有Picard工具集,所有Picard工具都能够使用GATK完成。
5、将过滤后的reads比对到参考基因组上,并将生成的sam文件转换成bam文件
## ID:S027为样品ID,SM:S027为样品名称,PL:illumina为测序平台,输入文件为参考基因组和质控过滤后的文件
bwa mem -t 10 -R '@RG\tID:S027\tSM:S027\tPL:illumina' genome.fasta S027.R1.fq.gz S027.R2.fq.gz | samtools sort -@ 4 -m 1G -o s027.sort.bam -
# bwa比对命令最好不要使用nohup,nohup的输出信息会被bwa输出到目标文件,可能会影响后续步骤。
# mem算法 -t 运行线程数 -R添加头部信息,-R一定不能省略或写错。
# bwa的mem的效率更高,且更加准确。
6、查看及统计比对结果
##查看bam文件
samtools view -h test.bam | less -S
##统计bam文件信息
samtools flagstat test.bam > test.bam.flagstat
7、去除重复
##去除bam文件中的重复
picard -Xmx4g MarkDuplicates I=test.sorted.bam O=test.sorted.rmdup.bam CREATE_INDEX=true REMOVE_DUPLICATES=true M=test.marked_dup_metrics.txt
# I 输入排序后的文件; O 输出文件; M 输出重复矩阵。
9、检测变异位点
##检测变异
nohup gatk --java-options "-Xmx10g" HaplotypeCaller -R genome.fasta -I test.sorted.rmdup.bam -ERC GVCF -O test.g.vcf &
10、合并变异位点
##生成合并列表
ls *.g.vcf > gvcf.list
##合并GVCF文件
nohup gatk --java-options "-Xmx10g" CombineGVCFs -R genome.fasta -V gvcf.list -O all.g.vcf &
11、提取所有样本变异位点
nohup gatk --java-options "-Xmx4g" GenotypeGVCFs -R genome.fasta -V all.g.vcf -O all.var.vcf &
12.1、分选SNP
## 从变异里提取所有SNP
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf --select-type SNP -O all.raw.snp.vcf
## 过滤不符合条件的SNP
gatk --java-options "-Xmx4g" VariantFiltration -R genome.fasta -V all.raw.snp.vcf --filter-expression "QD < 2.0 || MQ<40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name 'LowQualFilter' -O all.snp.fil.vcf
## 提取过滤后的SNP
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.snp.fil.vcf --exclude-filtered -O all.snp.filed.vcf
12.2、分选Indel
## 从变异里提取所有Indel
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf --select-type INDEL -O all.raw.indel.vcf
## 过滤不符合条件的Indel
gatk --java-options "-Xmx4g" VariantFiltration -R genome.fasta -V all.raw.indel.vcf --filter-expression "QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filter-name 'LowQualFilter' -O all.indel.fil.vcf
##提取过滤后的Indel
gatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.indel.fil.vcf --exclude-filtered -O all.indel.filed.vcf