本人通过无数个网上文章的指导，无数个萝卜的坑，终于摸索出了一整套完整的、傻瓜式的寻找全基因组的SNP实操流程。
仅以此文纪念我这俩个月的努力。
首先获得序列。序列的这一块我是已有的。所以也算不上训练。这只是一个流程化的东西。我做的也是蜜蜂的基因组。
在我所用的服务器里，samtools ， bwa ，gatk 都是预先构建好了全局变量。直接引用。如若没有设置好全局变量，请自行../../操作至软件目录下调用程序。
注意：我没有做BQSR的质控分析。
1、参考序列的索引的建立
samtools faidx reference.fasta （这是samtool所需的索引）
bwa index reference.fasta （这是bwa所需的索引）
java -jar picard.jar CreateSequenceDictionary REFERENCE=reference.fasta OUTPUT=reference.dict
加粗代表你的序列，斜体表示你的输出文件，点号的文件后缀不推荐更改。。
2、比对
bwa mem -t 16 -R '@RG\tID:X\tPL:illumina\tSM:X' refence.fasta R1.fasta R2.fasta | samtools view -Sb - > 1.bam
此处的 -t 16 代表的是处理线程，越多速度越快。，该步骤处理时间较长。似乎无法后台处理？自行判断。。
R1 R2代表测序得到的正反俩条序列，都要带上。
具体可以设置多少线程视服务器而定。可通过下列代码查询，设置相关值。
grep 'physical id' /proc/cpuinfo | sort -u 查看CPU个数
grep 'core id' /proc/cpuinfo | sort -u | wc -l 查看核心数量
grep 'core id' /proc/cpuinfo | sort -u | wc -l 查看线程
cat /proc/meminfo 查看服务器内存
3、排序
samtools sort -@ 16 -m 16G -O bam -o 1.sorted.bam 1.bam
将比对的bam结果进行排序 @后接的是线程， m后接的是内存
4、标记重复 （用时较长，可以选择后台跑）
gatk MarkDuplicates -I 1.sorted.bam -O 1.sorted.markdup.bam -M 1.sorted.markdup_metrics.txt
设置后台跑的只需在代码最前方加上 nohup 并在最后加上 & 就可设置后台跑了
5、创建对比索引文件
samtools index 1.sorted.markdup.bam
6、生成gvcf初步比对文件（处理时间长）
gatk HaplotypeCaller -R reference.fasta --emit-ref-confidence GVCF -I 1.sorted.markdup.bam -O 1.g.vcf
7、生成VCF文件
gatk GenotypeGVCFs -R reference.fasta -V 1.g.vcf -O test.vcf
到这里，SNP的初步筛选已经完成，后续就是VQSR的变异质控和过滤。
网友评论