参考教程：
https://cloud.tencent.com/developer/article/1054625
http://www.360doc.com/content/21/0714/12/76149697_986499282.shtml
http://www.bio-info-trainee.com/2163.html
一、统计基本覆盖信息
首先在linux里面运行：
samtools mpileup -f /data/zds209/database/cellranger/refdata-gex-GRCh38-2020-A/fasta/genome.fa /data/zds209/ssresult/bam/B8D.rm.bam| perl -alne '{$pos=int($F[1]/1000); $key="$F[0]\t$pos";$GC{$key}++ if $F[2]=~/[GC]/;$counts_sum{$key}+=$F[3];$number{$key}++;}END{print "$_\t$number{$_}\t$GC{$_}\t$counts_sum{$_}" foreach sort{$a<=>$b} keys %number}'
-f参数 输入有索引文件的fasta参考序列
用qsub提交以后，生成o文件就是具体内容。
image.png
下面这种也可以考虑使用，但是生成文件略有不同。
nohup time samtools mpileup P_jmzeng.final.bam |perl -alne '{$pos{$F[0]}++;$depth{$F[0]}+=$F[3]} END{print "$_\t$pos{$_}\t$depth{$_}" foreach sort keys %pos}' 1>coverage_depth.txt 2>coverage_depth.log &
关于samtool mpileup的详细说明参见教程：https://www.jianshu.com/p/a3791cf16474?ivk_sa=1024320u
https://cloud.tencent.com/developer/article/1441634
mpileup不使用-u或-g参数时，则不生成二进制的bcf文件，而生成一个文本文件(输出到标准输出)。该文本文件统计了参考序列中每个碱基位点的比对情况；该文件每一行代表了参考序列中某一个碱基位点的比对结果。上面命令就是把这个标准输出直接管道命令给perl 生成一个新的文件，或者输出到o文件。
二、做图可视化
把上面生成的o文件改名为coverage.txt，传到本地电脑。
在R里运行：
setwd("C:/temp/bed")
rm(list = ls())
file <- 'coverage.txt'
dat <- read.table(file, sep = "\t", fill=TRUE,stringsAsFactors = F)
colnames(dat)=c('chr','number','length','GC','counts')
keep_chr <- dat$chr %in% c(paste0('chr',c(1:22,'X','Y')))
dat <- dat[keep_chr,]
dat$depth <- dat$counts/dat$length
library(ggplot2)
head(dat)
# To change plot order of facet wrap,
# change the order of varible levels with factor()
dat$chr <- factor(dat$chr , levels =c(paste0('chr',c(1:22,'X','Y'))) )
png('coverage.png',height = 500,width = 2000)
p <- ggplot(dat,aes(number,depth))+geom_area(aes(fill=chr))+ylim(0, 10)
p <- p+facet_grid( ~ chr,scales = 'free_x',space = 'free_x')
print(p)
dev.off()
放个自己调整了的图，似乎比jimmy大神的看着好看一丢丢。
image.png
三、
这两个教程里面用了不同的分面。。（关于R的具体细节有待深入学习。。近期安排上）https://cloud.tencent.com/developer/article/1054625
用了 facet_wrap
https://cloud.tencent.com/developer/article/1928651
用了facet_grid
学习这两种函数的差别，教程如下：
https://www.jianshu.com/p/dd6bbbba117b
网友评论