对bam文件作基础统计
https://www.jianshu.com/p/4bc060bc6785
一只烟酒僧
0.078
2020.09.12 23:16:06
字数 76
阅读 245
参考链接https://www.jianshu.com/p/82ed6e27f571
一、有多少read比对到参考基因组上
使用软件：samtools flagstat或idxstats
参考samtools命令详解http://www.cnblogs.com/emanlee/p/4316581.html

命令：samtools flagstat mybam
结果解析：
3826122 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数
0 + 0 secondary
1658 + 0 supplementary
343028 + 0 duplicates
3824649 + 0 mapped (99.96% : N/A) #总体reads的匹配率
3824464 + 0 paired in sequencing #总共的reads数
1912442 + 0 read1 #reads1中的reads数
1912022 + 0 read2 #reads2中的reads数
3808606 + 0 properly paired (99.59% : N/A) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
3821518 + 0 with itself and mate mapped #paired reads中两条都比对到参考序列上的reads数
1473 + 0 singletons (0.04% : N/A) #单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
5882 + 0 with mate mapped to a different chr#paired reads中两条分别比对到两条不同的参考序列的reads数
4273 + 0 with mate mapped to a different chr (mapQ>=5) #paired reads中两条分别比对到两条不同的参考序列的reads数
命令:samtools idxstats mybam
结果：
chr1 248956422 53 0
chr2 242193529 30 0
chr3 198295559 30 0
chr4 190214555 26 0
chr5 181538259 14 0
chr6 170805979 18 0
chr7 159345973 18 0
chr8 145138636 22 0
chr9 138394717 20 0
chr10 133797422 30 0
chr11 135086622 30 0
chr12 133275309 18 0
chr13 114364328 14 0
chr14 107043718 18 0
chr15 101991189 21 0
chr16 90338345 18 0
chr17 83257441 22 0
chr18 80373285 8 0
chr19 58617616 26 0
chr20 64444167 10 0
chr21 46709983 2 0
chr22 50818468 16 0
chrX 156040895 24 0
chrY 57227415 0 0
chrM 16569 0 0
二、计算覆盖度coverage和测序深度depth
软件：
samtools depth
bedtools genomecov -ibam mybam -d >bedtools_genomecov_res

代码：bedtools genomecov -ibam mybam -d >bedtools_genomecov_res
结果：
chr1 1 0
chr1 2 0
chr1 3 0
chr1 4 0
chr1 5 0
chr1 6 0
chr1 7 0
chr1 8 0
chr1 9 0
chr1 10 0
chr1 11 0
chr1 12 0
代码：samtools depth mybam -a >samtools_depth_res
-a表示展示所有的pos，包括0深度的位点
结果：
chr1 1 0
chr1 2 0
chr1 3 0
chr1 4 0
chr1 5 0
chr1 6 0
chr1 7 0
chr1 8 0
chr1 9 0
chr1 10 0
chr1 11 0
chr1 12 0
chr1 13 0
chr1 14 0
chr1 15 0
chr1 16 0
chr1 17 0
chr1 18 0
chr1 19 0
chr1 20 0
chr1 21 0
chr1 22 0
chr1 23 0
chr1 24 0
chr1 25 0
chr1 26 0
chr1 27 0
chr1 28 0
chr1 29 0

我现在有四个样本的全基因组数据
我想统计如下数据：

1、计算1)全基因组平均测序深度 2）不同深度（1X 4X 10X 20X）下的测序覆盖度
for i in ls -d SRR*
do
{
bedtools genomecov -ibam $i/$i’.sorted.bam’|grep genome>$i/$i’.bedtools.genomecov.res’
}&
done
#接R的统计代码
info.list<-data.frame(matrix(0,nrow =length(grep(“SRR*”,list.files())) ,ncol =6 ))
n=0
for (i in list.files()[grep(“SRR*”,list.files())]) {
n=n+1
dir<-paste(i,"/",i,".bedtools.genomecov.res",sep = “”)
x<-read.table(dir,sep = “\t”)
x<-x[grep(“genome”,x$V1),]info.v<-vector()info.v[1]<-iinfo.v[2]<-sum(x$V3[-1]*x$V2[-1])/x$V4[1]
info.v[3]<-sum(x$V3[-1])/x$V4[1]
info.v[4]<-sum(x$V3[-c(1:4)])/x$V4[1]
info.v[5]<-sum(x$V3[-c(1:10)])/x$V4[1]
info.v[6]<-sum(x$V3[-c(1:20)])/x$V4[1]
info.list[n,]<-info.v
}
colnames(info.list)<-c(“ID”,“mean_depth”,“1X”,“4X”,“10X”,“20X”)

image.png
2、卡合适的阈值(质量值至少为30，深度至少为4，以及其它的硬阈值)的比对结果纳入高质量SNP范围（这里使用的gatk工具）

#!/bin/bash
for id in ls -d SRR*
do
{
input=$id/$id’_raw.vcf’
output=$id/$id’_snp.vcf’
gatk SelectVariants -select-type SNP -V $input -O $output
}&
done

for id in ls -d SRR*
do
{
input=$id/$id’_raw.vcf’
output=$id/$id’_indel.vcf’;
gatk SelectVariants -select-type INDEL -V $input -O $output
}&
done
wait

for id in ls -d SRR*
do
{
input=$id/$id’_snp.vcf’
output=$id/$id’_snp.filtered.vcf’
gatk VariantFiltration -V $input --filter-expression “DP<4 || QD < 2.0 || MQ < 40.0 || FS > 60.0 || SOR > 3.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” --filter-name “Filter” -O $output
}&
done 1>hard_filter_snp.log 2>hard_filter_snp.errlog

for id in ls -d SRR*
do
{
input=$id/$id’_indel.vcf’
output=$id/$id’_indel.filtered.vcf’
gatk VariantFiltration -V $input --filter-expression “DP<4 || QD < 2.0 || FS > 200.0 || SOR > 10.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0” --filter-name “Filter” -O $output
}&

done 1>hard_filter_indel.log 2>hard_filter_indel.errlog

wait

for id in ls -d SRR*
do
{
input1=$id/$id’_snp.filtered.vcf’
input2=$id/$id’_indel.filter.vcf’
output=$id/$id’.filtered.vcf’
gatk MergeVcfs -I $input1 -I $input2 -O $output
}&
done 1>merge_snp_indel.log 2>merge_snp_indel.errlog

#这里暂未考虑质量值的问题！！！

作者：一只烟酒僧
链接：https://www.jianshu.com/p/712e8aa30193
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

参考链接(注意这个是旧版，但是写的很详细！)：https://anjingwd.github.io/AnJingwd.github.io/2017/08/19/bedtools%E4%BD%BF%E7%94%A8%E6%95%99%E7%A8%8B%E8%AF%A6%E8%A7%A3/
参考链接：https://zhuanlan.zhihu.com/p/52322803

技巧一、使用bedtools bamtobed将bam文件转bed文件(基础)
参考：https://uteric.github.io/CHIPSEQ/bamtobed/

bedtools bamtobed -i SRR3022344.sorted.bam >test.bed
第一列：染色体位置
第二列：start
第三列：end
第四列：对应BAM文件的QNAME，包含测序平台，read name等信息
第五列：对应BAM文件的MAPQ，即比对质量
第六列：正负链
注意：
start1和start2起始坐标第一个碱基都为0，所以start=9, end=20表示碱基跨度是从第10位到第20位
染色体用.表示unknown；位置信息用-1表示unknown

image.png
技巧二、使用bedtools intersect计算两个或多个bed中的intersect区域(可接受多个文件类型bed/gff/vcf/bam)
bedtools intersect [OPTIONS] -a <bed/gff/vcf/bam> -b <bed/gff/vcf/bam>

-wa参数可以报告出原始的在A文件中的feature
-wb参数可以报告出原始的在B文件中的feature
-c参数可以报告出两个文件中的overlap的feature的数量
-wo 返回overlap碱基数
-v 返回非overlap区间
-s 相同链上的feature

#例题请看参考链接！
注意，自己生成测试bed文件，都必须用tab键分割，否则会报错！！

案例一：包含着染色体位置的两个文件，分别记为A文件和B文件。分别来自于不同文件的染色体位置的交集是什么？

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed

chr1 15 20

案例二：包含着染色体位置的两个文件，分别记为A文件和B文件。求A文件中哪些染色体位置是与文件B中的染色体位置有overlap.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -wa

chr1 10 20

案例三：包含着染色体位置的两个文件，分别记为A文件和B文件。求A文件中染色体位置与文件B中染色体位置的交集，以及对应的文件B中的染色体位置.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -wb

chr1 15 20 chr1 15 25

案例四（经用）： 包含着染色体位置的两个文件，分别记为A文件和B文件。求对于A文件的染色体位置是否与文件B中的染色体位置有交集。如果有交集，分别输入A文件的染色体位置和B文件的染色体位置；如果没有交集，输入A文件的染色体位置并以’. -1 -1’补齐文件。

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 25

$ bedtools intersect -a A.bed -b B.bed -loj

chr1 10 20 chr1 15 25

chr1 30 40 . -1 -1

案例五： 包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置，以及overlap的长度.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -wo

chr1 10 20 chr1 15 20 5

chr1 10 20 chr1 18 25 2

案例六： 包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，如果和B文件中染色体位置有overlap,则输出在A文件中染色体位置和在B文件中染色体位置，以及overlap的长度；如果和B文件中染色体位置都没有overlap,则用’. -1-1’补齐文件

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -wao

chr1 10 20 chr1 15 20 5

chr1 10 20 chr1 18 25 2

chr1 30 40 . -1 -1

案例七： 包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，输出在A文件中染色体位置和有多少B文件染色体位置与之有overlap.

$ cat A.bed

chr1 10 20

chr1 30 40

$ cat B.bed

chr1 15 20

chr1 18 25

$ bedtools intersect -a A.bed -b B.bed -c

chr1 10 20 2

chr1 30 40 0

案例八(常用)： 包含着染色体位置的两个文件，分别记为A文件和B文件。对于A文件中染色体位置，输出在A文件中染色体位置和与B文件染色体位置至少有X%的overlap的记录。

$ cat A.bed

chr1 100 200

$ cat B.bed

chr1 130 201

chr1 180 220

$ bedtools intersect -a A.bed -b B.bed -f 0.50 -wa -wb

chr1 100 200 chr1 130 201

#一个2G的bam文件约需要15-20分钟！！！
注意：当用bedtools intersect 处理大文件时比较耗内存，有效的方法是对A和B文件按照染色体名字(chromosome)和位置(position)排序，重新intersect

bedtools sort [OPTIONS] -i <bed/gff/vcf>
bedtools sort -chrThenSizeA -i test.bed
-sizeA Sort by feature size in ascending order.
-sizeD Sort by feature size in descending order.
-chrThenSizeA Sort by chrom (asc), then feature size (asc).
-chrThenSizeD Sort by chrom (asc), then feature size (desc).
-chrThenScoreA Sort by chrom (asc), then score (asc).
-chrThenScoreD Sort by chrom (asc), then score (desc).
技巧三：使用bedtools genomecov染色体和全基因组覆盖度计算（可以用来做深度计算）
单个输入bed文件（-i指定）和genome files；如果输入为bam(-ibam指定)文件，则不需要genome files

bedtools genomecov [OPTIONS] -i <bed/gff/vcf> -g
-ibam The input file is in BAM format
Note: BAM must be sorted by position

示例
$ cat ranges-cov-sorted.bed
chr1 4 9
chr1 1 6
chr1 8 19
chr1 25 30
chr2 0 20

$ cat cov.txt （染色体及每条染色体总碱基数）
chr1 30
chr2 20

bedtools genomecov -i ranges-cov-sorted.bed -g cov.txt
chr1 0 7 30 0.233333 1
chr1 1 20 30 0.666667
chr1 2 3 30 0.1
chr2 1 20 20 1 2
genome 0 7 50 0.14 3
genome 1 40 50 0.8
genome 2 3 50 0.06
#name 覆盖次数 覆盖碱基数 总碱基数 覆盖度
#同时计算单染色体和全基因组覆盖度

如果输入的是-ibam bam 文件，则输出结果为
名称 覆盖次数 该次数下的碱基数 该染色体长度 染色体覆盖度
chr1 0 248951704 248956422 0.999981
chr1 1 4136 248956422 1.66134e-05
chr1 2 582 248956422 2.33776e-06
chr2 0 242190850 242193529 0.999989
chr2 1 2358 242193529 9.73602e-06
chr2 2 321 242193529 1.32539e-06
chr3 0 198292860 198295559 0.999986
chr3 1 2398 198295559 1.20931e-05
chr3 2 301 198295559 1.51794e-06
chr4 0 190212444 190214555 0.999989
chr4 1 1622 190214555 8.52721e-06
chr4 2 489 190214555 2.57078e-06
chr5 0 181536919 181538259 0.999993
chr5 1 1280 181538259 7.05086e-06
chr5 2 60 181538259 3.30509e-07
chr6 0 170804532 170805979 0.999991
chr6 1 1095 170805979 6.41078e-06
chr6 2 352 170805979 2.06082e-06
chr7 0 159344307 159345973 0.99999

chr1 248951704+4136+582=248956422

2G的bam文件大概的时间消耗
real 3m29.649s
user 3m19.071s
sys 0m9.387s
genomecov也会对全基因组的按照不同深度做统计！！（有用啊）
-u Write the original A entry once if any overlaps found in B
注意-u参数，1、使用后相当于进入-wa模式，2、若A与B有相同重复，则会去重
image.png
技巧四、使用bedtools coverage计算染色体给定区间的深度和覆盖度，输入文件可以是bam

bedtools coverage [OPTIONS] -a <bed/gff/vcf> -b <bed/gff/vcf>
其中-a指定interval文件，即你想看的染色体区间 -b指定的是你比对的结果或bed等文件

示例
$ cat A.bed
chr1 0 100
chr1 100 200
chr2 0 100

$ cat B.bed
chr1 10 20
chr1 20 30
chr1 30 40
chr1 100 200

$ bedtools coverage -a A.bed -b B.bed
A的名称 起始 终止 B在A中匹配次数 匹配的总长度 该区域总长度 比例
chr1 0 100 3 30 100 0.3000000
chr1 100 200 1 100 100 1.0000000
chr2 0 100 0 0 100 0.0000000
结果解释：前三列是interval文件的信息，后四列为统计信息

作者：一只烟酒僧
链接：https://www.jianshu.com/p/61bdad6288b3
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。