「来绘图吧!」链特异性测序的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.

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,529评论 5 475
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,015评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,409评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,385评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,387评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,466评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,880评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,528评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,727评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,528评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,602评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,302评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,873评论 3 306
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,890评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,132评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,777评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,310评论 2 342

推荐阅读更多精彩内容