以釀酒酵母為例，用到的數(shù)據(jù)如下，
隨便找一個鏈特異性測序的fastq用常見比對軟件跑一下都行。
有一些數(shù)據(jù)我放在github，可參考：https://github.com/Youpu-Chen/Myscripts/tree/miniproject/mini-project/strand-specific-RNAseq
- Saccharomyces cerevisiae的基因組文件
- Saccharomyces cerevisiae的鏈特異性測序數(shù)據(jù)
上游處理
稍微解釋一下為什么要這樣做。
普通的RNA-Seq建庫測序方式，反轉(zhuǎn)錄得到的cDNA，連接上的adaper不具有特異性，只是為了將目標序列連接到oligo上。最終的reads alignment無法區(qū)分開比對到某一個區(qū)域（包括positive strand和negative strand的一段基因組序列）的reads，來自該區(qū)域正負鏈的真實性。
Note：真實性，指這段區(qū)域的正負鏈真的在轉(zhuǎn)錄過程中有表達reads，而不是由于反轉(zhuǎn)錄建庫而引起的。
而鏈特異性測序為例（以加入dUTP的方法為例，也是使用最廣泛的一種方法）。shortgun得到目標mRNA序列之后，第二輪直接加入dUTP作為標記，之后再把含有這樣標記的序列給拿掉。那么上述的操作，就保證了最終基于read alignments的定量結(jié)果，能夠精確到正負鏈。
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即可。比如與人類免疫有關(guān)的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)
}
結(jié)果圖如下，
寫在后面
最近放開，陽了好多人，但是我還是堅持能不陽就不陽的原則。結(jié)果跑來跑去的身體搞垮的，陽倒是沒陽，人感冒了。但是剛好混到一點顆粒和靠曾老師給的一箱橙子茍活。
2022年很艱難，但是只要選對方向，不斷努力就好！
年末了，好想寫一篇長長的總結(jié)，來感謝我今年遇到的人和事。
這個其實是我期末大作業(yè)的一部分，若有同學(xué)有幸看到了，歡迎和我討論。
還有就是關(guān)于RNA-Seq、strand-specific RNA-Seq對最后的DEG分析有沒有影響，統(tǒng)計檢驗的power怎么樣，感興趣的同學(xué)可以拿負二項分布和泊松分布擬合試試，我懶，我不寫了。
參考資料
[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.
