用组装出来的线粒体序列call SNP

想用线粒体基因组建一个树，看一下系统发生关系。

一个思路是从文献中看到的，将线粒体序列的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的图