「来绘图吧！」链特异性测序的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)
}
结果图如下，
写在后面
最近放开，阳了好多人，但是我还是坚持能不阳就不阳的原则。结果跑来跑去的身体搞垮的，阳倒是没阳，人感冒了。但是刚好混到一点颗粒和靠曾老师给的一箱橙子苟活。
2022年很艰难，但是只要选对方向，不断努力就好！
年末了，好想写一篇长长的总结，来感谢我今年遇到的人和事。
这个其实是我期末大作业的一部分，若有同学有幸看到了，欢迎和我讨论。
还有就是关于RNA-Seq、strand-specific RNA-Seq对最后的DEG分析有没有影响，统计检验的power怎么样，感兴趣的同学可以拿负二项分布和泊松分布拟合试试，我懒，我不写了。
参考资料
[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.
