bam格式文件处理大全（五）

sam文件是短序列比对生成的文件，是二代测序中最核心的文件。在RNAseq，变异检测等分析中，都需要首先生成sam文件格式。bam文件是sam格式的二进制格式，转换为二进制之后，可以减小文件的存储。掌握sam/bam文件的操作是处理二代测序数据的非常重要的内容，例如sam与bam的转换，排序，建立索引，reads计数等等操作。

21 mpileup

mpileup是将reads堆叠起来，作用是展示每一个位点比对的细节，然后可以通过这个细节来推测比对情况，pileup是之前maq工具查找SNP重要的步骤，目前已经逐渐将mpileup功能转移到bcftools中了。

 samtools mpileup --reference ref/ref.fna   A1.sorted.bam 

22 MarkDuplication

Dupliacation reads会对变异检测造成干扰，得到一些假阳性的结果，因此，需要将这些reads去除掉，只留一份即可。可以在比对之前去除，但是比较消耗内存。也可以在比对之后进行标记。这一步骤只是在每一行比对结尾出添加一些CIGAR标志，并不过滤数据。samtools可以标记Duplication，也可以去除掉reads，GATK也可以进行标记。

gatk MarkDuplicates -I A1.sorted.bam -M Sample1.markdup_metrics.txt -O A1.sorted.markdup.bam  
samtools index A1.sorted.markdup.bam  

23 BQSR

BQSR是Base Quality Score Recalibratio的简称，是利用GATK做人全基因组分析的必要步骤。首先利用BaseRecalibrator进行机器学习，然后利用ApplyBQSR应用于比对结果。但是需要注意，BQSR至少要达到100万个突变位点，因此，外显子和目标区域捕获一般达不到这个要求，也就做不了这步奏。

#进行学习
 gatk BaseRecalibrator \
         -R Homo_sapiens_assembly38.fasta \
         -I A1.sorted.markdup.bam \
         --known-sites 1000G_phase1.snps.high_confidence.hg38.vcf.gz \
         --known-sites Mills_and_1000G_gold_standard.indels.hg38.vcf.gz \
         --known-sites dbsnp_146.hg38.vcf.gz \
         -O A1.sorted.markdup.recal_data.table >bqsr.log
#进行BQSR
gatk ApplyBQSR \
         --bqsr-recal-file A1.sorted.markdup.recal_data.table \
         -R Homo_sapiens_assembly38.fasta \
         -I A1.sorted.markdup.bam \
         -O A1.sorted.markdup.BQSR.bam
#对新的bam创建索引
samtools index A1.sorted.markdup.BQSR.bam

24 利用bcftools进行SNP检测

可以使用bcftools直接来筛选SNP，输入排序并建立索引的bam即可，如果能做Mark Duplication则更好了。

#利用bcftools进行变异检测     
bcftools mpileup -f ref.fna A1.sorted.bam | bcftools call -c -v -o A1.bcftool.vcf  

25 利用freebayes进行snp检测

与bcftools类似，使用freebayes来筛选SNP也非常容易，输入排序并建立索引的bam。

freebayes  -f ref.fna -C 5 A1.sorted.bam -v A1.freebayes.vcf 

---------- END ----------

（添加作者微信，请注明单位姓名）

您可能还会感兴趣的