想用线粒体基因组建一个树，看一下系统发生关系。
一个思路是从文献中看到的，将线粒体序列的13个蛋白编码基因分别先mafft比对，然后再串联到一起建树
另一个思路用线粒体基因组call SNP，详情见下图
第一种做法是用个体的重测序文件比对到参考(线粒体)基因组上；
第二种也是本文所用的，是用重测序文件组装出来的线粒体序列比对到参考(线粒体)基因组上，但是这种方法靠性存疑。
0.用getorganelle组装线粒体序列
详情见https://github.com/Kinggerm/GetOrganelle
1.比对到参考基因组上(bwa mem，用bwa aln的话call不出来，原因不明)
for i in {xxx} 这部分下面省略
do
bwa mem -t 4 -R "@RG\tID:$i\tPL:illumina\tSM:$i\tLB:$i" ref.fasta $i.fasta | samtools view -bS- > $i.bam
done
tID和tSM和tLB要保持一致，这里设置的是 个体名，否则call不出来，原因暂不明。tPL即测序平台。
2.排序、标记重复序列、建索引
samtools sort -@ 3 $i.bam -o $i.sorted.bam
gatk MarkDuplicates -I $i.sorted.bam -O $i.sorted.markdup.bam -M $i.sorted.markdup_metrics.txt
samtools index $i.sorted.markdup.bam
http://www.360doc.com/content/19/1224/14/68068867_881793271.shtml 标记重复序列的意义与作用，-M参数即把标记到的重复序列输出
3.call SNP
samtools mpileup 1.bam 2.bam… -f ref.fasta –gD –o test.bcf
bcftools view test.bcf > test.vcf 这两行代码的作用和 bcftools mpileup等效
bcftools call test.vcf -c -v -o end.vcf
这行代码的作用见下图，参考自https://www.bioinfo-scrounger.com/archives/248/
另外，用gatk就call不出来，所以换成了bcftools
最后得到了end.vcf便是结果，因为这种做法用到的线粒体序列数据很小，一个个体一百多kb，所以没有进行过滤
画了一个PCA的图
网友评论