根据bam文件按特定区间长度统计测序深度并画折线图展示结果

第一步使用bwa将测序结果比对到参考基因组，获得sam格式文件

bwa index pra2.fasta
bwa mem -t 8 -R '@RG\tID:foo\tSM:sample1' pra2.fasta filtered_1_paired.fq filtered_2_paired.fq -o pra.sam

sam转换为bam并排序

samtools view -S -b -o pra.bam pra.sam
samtools sort pra.bam -o pra.sorted.bam

统计测序深度用到的软件是BamDeal

软件对应的github主页 GitHub - BGI-shenzhen/BamDeal: BamDeal: a comprehensive toolkit for bam manipulation

软件下载链接 Download (github.com)

软甲解压出来后在bin目录下有一个可执行程序，可以直接使用，按照GitHub主页的内容首先赋予软件可执行程序

首先是解压

tar -xzvf BamDeal-0.24.tar.gz

切换至bin目录下

cd BamDeal-0.24/bin/

赋予权限

chmod 755 BamDeal_Linux

BamDeal_Linux 这个文件会由原来的的白色变为绿色

统计测序深度

BamDeal_Linux statistics Coverage -i pra.sorted.bam -r pra2.fasta -o depthpra -w 1000

-i 参数指定bam文件
-r 参数指定参考序列
-o 参数指定输出文件的前缀
-w 参数指定区间长度，默认好像是10000

测序深度的结果文件是depthpra.DepthGC.gz 同时里面还有一个数据是GC含量
这个文件的内容

##Bigin Depth(X)    GC(%)
#>Chr1  WindowsSite:    1000
1   408.26  35.10
1001    425.90  36.30
2001    454.05  32.20
3001    450.86  31.30
4001    407.55  27.20
5001    365.52  29.50
6001    401.80  29.00
7001    397.63  29.80
8001    405.23  30.40
9001    404.18  29.00
10001   351.01  27.60
11001   458.65  40.80
12001   449.08  34.30
13001   407.47  35.20
14001   414.81  28.60
15001   461.80  35.40
16001   434.95  35.40
17001   440.03  38.10
18001   444.49  32.60
19001   476.05  34.90
20001   464.12  40.70
21001   465.50  38.30
22001   466.07  40.40
23001   494.14  35.50
24001   501.55  38.30
25001   452.97  41.00
26001   502.83  37.90
27001   510.33  30.40
28001   435.41  32.60
29001   449.89  27.60
30001   472.28  29.50
31001   449.31  30.90
32001   427.53  38.00
33001   440.45  27.10
34001   423.36  35.20
35001   433.87  44.30
36001   428.36  43.70
37001   429.49  36.10
38001   400.11  35.10
39001   496.11  38.50
40001   488.33  41.50
41001   433.54  42.70
42001   467.28  42.90
43001   434.32  40.40
44001   421.64  33.40
45001   447.50  36.90
46001   426.81  31.60
47001   419.58  36.90
48001   281.95  23.60
49001   456.64  33.50
50001   403.85  29.70
51001   441.83  36.30
52001   429.97  34.10
53001   287.46  26.10
54001   469.80  35.30
55001   557.67  42.50
56001   493.17  35.10
57001   465.26  40.50
58001   439.33  37.20
59001   402.75  29.90
60001   405.90  35.50
61001   363.55  27.20
62001   424.82  33.80
63001   458.71  32.00
64001   447.40  40.30
65001   449.99  29.20
66001   460.24  37.30
#>Chr2  WindowsSite:    1000
1   473.05  33.40
1001    373.05  30.50
2001    411.34  33.20
3001    465.39  33.90
4001    432.12  31.60
5001    434.39  33.30
6001    438.59  37.90
7001    385.72  43.80
8001    476.58  35.70
9001    454.43  33.50
10001   493.96  37.90
11001   463.05  36.50
12001   518.58  31.20
13001   477.54  37.20
14001   477.00  31.70
15001   443.77  37.10
16001   478.95  30.70
17001   520.21  33.70
18001   258.27  37.60
19001   0.00    38.90
20001   0.00    37.20
21001   0.00    38.10
22001   0.00    37.30
23001   0.00    34.10
24001   0.00    36.10
25001   0.00    37.60
#>Chr3  WindowsSite:    1000
1   0.00    41.70
1001    0.00    38.90
2001    0.00    39.20
3001    0.00    37.70
4001    0.00    39.00
5001    0.00    37.70
6001    0.00    40.10
7001    0.00    38.10
8001    0.00    45.80
9001    0.00    55.80
10001   0.00    53.60
11001   0.00    48.90
12001   0.00    49.70
13001   0.00    55.20
14001   0.00    56.10
15001   0.00    52.50
16001   0.00    45.30
17001   0.92    36.70
18001   352.57  31.10
19001   520.56  34.40
20001   464.88  28.40
21001   356.30  24.00
22001   349.78  21.70
23001   482.00  33.90
24001   473.14  32.30
25001   485.40  36.00
26001   471.24  33.20
27001   488.30  29.90
28001   503.08  31.10
29001   445.33  30.70
30001   466.26  36.10
31001   501.01  38.10
32001   471.55  29.60
33001   498.57  28.70
34001   484.16  26.70
35001   467.88  24.90
36001   404.52  29.90
37001   329.56  32.10
38001   0.02    37.10
39001   0.00    45.30
40001   0.00    53.10
41001   0.00    56.60
42001   0.00    54.60
43001   0.00    49.80
44001   0.00    48.10
45001   0.00    54.60
46001   0.00    55.50
47001   0.00    45.20
48001   0.00    37.90
49001   0.00    40.10
50001   0.00    37.50
51001   0.00    39.70
52001   0.00    37.30
53001   0.00    39.70
54001   0.00    38.40
55001   273.93  42.10
56001   206.32  37.50
57001   0.00    36.30
58001   0.00    34.30
59001   0.00    35.80
60001   0.00    38.90
61001   0.00    37.10
62001   0.00    39.50

这个数据导出我们自己画图，这个数据怎们整理我暂时还想不到好的办法，分成3个文件分别读入吧

df1<-read.table("chr1.txt",header=T,sep="\t")
df2<-read.table("chr2.txt",header=T,sep="\t")
df3<-read.table("chr3.txt",header=T,sep="\t")
df1$group<-"Chr1"
df1$x<-1:dim(df1)[1]
colnames(df1)<-c("Chr","y","gc","group","x")

df2$group<-"Chr2"
df2$x<-1:dim(df2)[1]
colnames(df2)<-c("Chr","y","gc","group","x")

df3$group<-"Chr3"
df3$x<-1:dim(df3)[1]
colnames(df3)<-c("Chr","y","gc","group","x")

df<-rbind(df1,df2,df3)
df
library(ggplot2)
ggplot(df,aes(x=x,y=y))+
  geom_point(aes(color=group))+
  geom_line(aes(color=group))+
  facet_wrap(~group,ncol=1,nrow=3,
             scales = "free")+
  theme_bw()+
  theme(legend.position = "none")

image.png

欢迎大家关注我的公众号
小明的数据分析笔记本