本人通过无数个网上文章的指导,无数个萝卜的坑,终于摸索出了一整套完整的、傻瓜式的寻找全基因组的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的变异质控和过滤。