「来绘图吧！」链特异性测序的read depth绘制

以酿酒酵母为例，用到的数据如下，
随便找一个链特异性测序的fastq用常见比对软件跑一下都行。
有一些数据我放在github，可参考：https://github.com/Youpu-Chen/Myscripts/tree/miniproject/mini-project/strand-specific-RNAseq

• Saccharomyces cerevisiae的基因组文件
• Saccharomyces cerevisiae的链特异性测序数据

上游处理

稍微解释一下为什么要这样做。
普通的RNA-Seq建库测序方式，反转录得到的cDNA，连接上的adaper不具有特异性，只是为了将目标序列连接到oligo上。最终的reads alignment无法区分开比对到某一个区域（包括positive strand和negative strand的一段基因组序列）的reads，来自该区域正负链的真实性。

Note：真实性，指这段区域的正负链真的在转录过程中有表达reads，而不是由于反转录建库而引起的。

而链特异性测序为例（以加入dUTP的方法为例，也是使用最广泛的一种方法）。shortgun得到目标mRNA序列之后，第二轮直接加入dUTP作为标记，之后再把含有这样标记的序列给拿掉。那么上述的操作，就保证了最终基于read alignments的定量结果，能够精确到正负链。

Note：普通的建库测序方法，就相当于放大镜，直接将大家的表达量都乘以2，最后再进行different gene expression。

上游分析代码如下，

# make chromosome bed file
samtools faidx genome.fa
awk '{print $1"\t"$2}' genome.fa.fai > genome.txt
bedtools makewindows -g genome.txt -w 1000 > windows.bed

# generate depth
samtools view -f 16 -b wtssRNA_seq.fastp.sorted.sam > wtssRNA_seq.fastp.sorted.rev.bam
samtools view -F 16 -b wtssRNA_seq.fastp.sorted.sam > wtssRNA_seq.fastp.sorted.fwd.bam
bedtools coverage -a windows.bed -b input.sorted.rev.bam > rev.depth.txt
bedtools coverage -a windows.bed -b input.sorted.fwd.bam > fwd.depth.txt

下游画图

rm(list=ls())

# Load datasets
fwd.df <- read.delim('fwd.depth.txt', sep = '\t', header = FALSE)
names(fwd.df) <- c('chr', 'start', 'end', 'read.counts', 'base.number', 'window.size', 'average.coverage')
head(fwd.df)

rev.df <- read.delim('rev.depth.txt', sep = '\t', header = FALSE)
names(rev.df) <- c('chr', 'start', 'end', 'read.counts', 'base.number', 'window.size', 'average.coverage')
head(rev.df)


# Load packages
library(ggplot2)
library(ggpubr)


# plot
if(F){
  ### chrIV，选择感兴趣的chromosome即可。比如与人类免疫有关的chrVI就是一个非常好的可视化对象
  fwd.chrI <- fwd.df[which(fwd.df$chr == "chrIV"), ]
  rev.chrI <- rev.df[which(rev.df$chr == "chrIV"), ]
  p1 <- ggplot(data = fwd.chrI, aes(x=(start+end)/2, y = read.counts)) + 
    # geom_line() + 
    geom_area(colour = "white", fill="red", alpha=0.5, position="identity") +
    scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0), 
                                                              breaks = c(0, 250, 500, 750), 
                                                              limits = c(0, 750)) + 
    theme_classic() + 
    xlab("Chromosome IV (Bp)") +
    # ylab("Read Counts/Per 1kb") + 
    # xlab("") + 
    ylab("") +
    # theme(
    #   axis.title.x = element_text(vjust=2)
    # ) + 
    annotate("text", x = 1350000, y = 600, label = "Positive Strand", size=5)
  p1
  
  
  p2 <- ggplot(data = rev.chrI, aes(x=(start+end)/2, y = read.counts)) + 
    # geom_line() + 
    geom_area(colour = "white", fill="blue", alpha=0.5, position="identity") + 
    theme_classic() + 
    scale_x_continuous(expand = c(0, 0)) + scale_y_continuous(expand = c(0, 0), 
                                                              breaks = c(0, 250, 500, 750)) +  
    scale_y_reverse(breaks = c(0, 250, 500, 750),
                    limits = c(750, 0)) + 
    geom_hline(yintercept = 0) + 
    xlab("") + 
    ylab("") +
    annotate("text", x = 1350000, y = 600, label = "Negative Strand", size=5)
  p2
  
  
  figure <- ggarrange(p1, p2,
                      ncol = 1, nrow = 2)
  figure
  
  
  anno.fig <- annotate_figure(
    figure,
    left = text_grob("Read Depth (Per 1kb)",
                     color = "black", rot = 90)
  )
  anno.fig
  
  #save 
  ggsave('chrIV-readdepth.pdf', anno.fig, height = 6, width = 10)
}

结果图如下，

参考资料

[1] Luan, H., Meng, N., Fu, J., Chen, X., Xu, X., Feng, Q., Jiang, H., Dai, J., Yuan, X., Lu, Y. and Roberts, A.A., 2014. Genome-wide transcriptome and antioxidant analyses on gamma-irradiated phases of Deinococcus radiodurans R1. PloS one, 9(1), p.e85649.