1、测序质量评估
## conda安装fastqcfastqc *.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.fastabwa 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.fastaGATK(全称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
## 从变异里提取所有SNPgatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf --select-type SNP -O all.raw.snp.vcf## 过滤不符合条件的SNPgatk --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 ## 提取过滤后的SNPgatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.snp.fil.vcf --exclude-filtered -O all.snp.filed.vcf
12.2、分选Indel
## 从变异里提取所有Indelgatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.var.vcf --select-type INDEL -O all.raw.indel.vcf## 过滤不符合条件的Indelgatk --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 ##提取过滤后的Indelgatk --java-options "-Xmx4g" SelectVariants -R genome.fasta -V all.indel.fil.vcf --exclude-filtered -O all.indel.filed.vcf
如果觉得《GWAS数据分析流程—SNP Indel提取》对你有帮助，请点赞、收藏，并留下你的观点哦！
