生物数据格式 - SAM/BAM

格式

sam文件是短序列比对生成的文件，是二代测序中最核心的文件。在RNAseq，变异检测等分析中，都需要首先生成sam文件格式。bam文件是sam格式的二进制格式，转换为二进制之后，可以减小文件的存储。sam文件的序列格式如下： sam/bam示例.png

在header之后，每一行是一条record，也就是一条read比对的结果，信息包括标配的十一列和可能有的第十二列：
第一列：reads ID

第二列：是flag标识符，是一些代表比对模式的二进制数字之和，数字的具体含义如下： flag含义.png

第三列：比对到参考序列上的染色体号。
第四列：在参考序列上的位置
第五列：比对的质量值，MAPQ
第六列：代表比对结果的CIGAR字符串
第七列：比对到的染色体号，若是没有比对到，则是*
第八列：比对到参考序列上的第一个碱基位置
第九列：Template的长度
第十列：为read的序列
第十一列：为ASCII码格式的序列质量

第十二列：元信息是第六列的详细信息，展示了比对的具体细节，一般包括了RG信息、missmatch信息、二次比对信息等。不过这些标识符是供下游软件处理，并不需要人工来解读，具体含义如下： image.png

获取

许多比对软件如bwa，bowie2，tophat2，subread，minimap2等等都可以生成sam格式。软件将fastq格式的测序数据比对回fasta格式的参考序列上，就得到了储存了比对信息的sam文件。

# bwa建立索引
bwa index -a is ref.fna
#bwa mem比对
bwa mem -t 4 -R '@RG\tID:A1\tPL:illumina\tSM:A1'  ref.fa  A1_1.fastq.gz A1_2.fastq.gz >A1.sam

bwa默认生成的是sam格式，而绝大部分的软件需要bam格式，所以需要将sam转换为二进制的bam，可以节省大量存储。

samtools view -O bam -o A1.bam A1.sam  #sam转bam
samtools view A1.bam | less -S #查看bam

有些文件格式不正确或者内容不完整的bam文件运行samtools会报错，这个时候就可以使用samtools自带的quickcheck功能进行验证。

samtools quickcheck *.bam  && echo 'all ok' || echo 'fail!'

处理

文件与序列处理

##排序
samtools sort -@ 4 -m 12G -O bam -o A1.sorted.bam A1.sam  
##建立索引
samtools index A1.sorted.bam  
##统计
samtools stats  A1.sorted.bam  > A1.stats #比对结果统计
plot-bamstats -p test  A1.stats  #并绘图
samtools flags 83 #查看flag含义
samtools flagstat A1.sorted.bam #flag统计
samtools idxstats male.sorted.bam | head #单条染色体比对统计
samtools bedcov cancer_panel.bed A1.sorted.bam #目标区域(bed)统计
samtools depth A1.sorted.bam#深度统计
##过滤
samtools view -f 4 A1.sorted.bam #将没有比对上的reads筛选出来
samtools veiw -F 4 A1.sorted.bam #将比对上的reads输出出来
##输出比对fastq与fasta
samtools fastq A1.sorted.bam -1 A.1.fq.gz -2 A.2.fq.gz -c 6 #输出比对fastq
samtools fasta A1.sorted.bam #输出fasta
##tview查看每个位点细节
samtools tview A1.sorted.bam #排序建立索引的reads
samtools tview A1.sorted.bam ref.fna #排序建立索引的参考序列
##可视化
      #tablet
      #IGV
##mpileup
samtools mpileup --reference ref/ref.fna   A1.sorted.bam #pileup是之前maq工具查找SNP重要的步骤，目前已经逐渐将mpileup功能转移到bcftools中

功能相关处理

##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 #GATK进行标记
samtools index A1.sorted.markdup.bam #samtools进行标记
##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
##SNP检测   
bcftools mpileup -f ref.fna A1.sorted.bam | bcftools call -c -v -o A1.bcftool.vcf #bcftools进行检测  
freebayes  -f ref.fna -C 5 A1.sorted.bam -v A1.freebayes.vcf #freebayes进行检测
### GATK进行检测，输入数据为排序建立索引的bam，如果是人全基因组，则为排序，标记Duplication并经过BQSR的bam，并建立索引。
gatk HaplotypeCaller --emit-ref-confidence GVCF -R ref.fna -I A1.sorted.bam -O  A1.g.vcf.gz #HaplotypeCaller：Call germline SNPs and indels via local re-assembly of haplotypes 
gatk GenotypeGVCFs -R ref.fna A1.g.vcf.gz -O A1.HC.vcf.gz #GenotypeGVCFs：Perform joint genotyping on gVCF files produced by HaplotypeCaller
##SV检测
delly call -g ref.fna -o A11.delly.sv.bcf -n A11.sorted.bam  #delly进行检测
delly filter -f germline -p -q 20 A1.delly.sv.bcf -o A1.delly.sv.filter.bcf #germline突变过滤
————————————————————————————————————————————————————————————————————————————————————————————
samtools view -b -F 1294A1.sorted.bam  | samtools sort - > A1.discordants.sorted.bam #1 挑选flag值为1294(非正常范围内的比对)的reads 
samtools view -h A1.sorted.bam | extractSplitReads_BwaMem -i stdin | samtools view -Sb - | samtools sort -> A.splitters.sorted.bam  #2 挑选split比对的reads
lumpyexpress -BA1.sorted.bam -S A1.discordants.sorted.bam -D A1.splitters.sorted.bam -o A1.lumpu.sv.vcf  #3 利用lumpy找SV
##CNVnator进行CNV检测
cnvnator -root A1.root -tree A1.sorted.bam -unique #1.提取mapping信息  
cnvnator -root A1.root -his 100  -d ref.fna #2.生成质量分布图HISTOGRAM   
cnvnator -root A1.root -stat 100 #3.生成统计结果  
cnvnator -root A1.root -partition 100  #4.RD信息分割partipition   
cnvnator -root A1.root -call 100 > A1.cnvnator.vcf #5.变异检出

———以上纯属个人理解与记录