GRO-seq 分析
第一篇就以最近做過(guò)的GRO-seq分析為例來(lái)展開(kāi)吧。
GRO-seq的原理和用處
GRO-seq 全稱 Global run-on sequencing, 是一種新生轉(zhuǎn)錄本RNA的測(cè)序方法（Core et al., 2008）。它能夠map基因組上RNA polymerase正在轉(zhuǎn)錄的位置，并且定量。在RNA 合成時(shí)，摻入NTP類似物bromouridine（BrU）, 新合成的RNA 便帶上了BrU，BrU可以通過(guò)抗體被capture下來(lái)。RNA反轉(zhuǎn)錄成cDNA后，可以被建庫(kù)測(cè)序，從而得知RNA合成的位置和定量。這樣能夠準(zhǔn)確地反映出某一時(shí)刻RNA polymerase的位置和活動(dòng)。可不要小看這個(gè)方法，通過(guò)這個(gè)方法，發(fā)現(xiàn)了在人類基因中，RNA polymerase 停在promoter附近的現(xiàn)象；還發(fā)現(xiàn)了RNA polymerase 在終止轉(zhuǎn)錄時(shí)，會(huì)越過(guò)基因的3‘end下游（pre-mRNA 3' end cleavage site）。另外，也發(fā)現(xiàn)enhancer也可以產(chǎn)生短的不穩(wěn)定RNA。關(guān)于這些發(fā)現(xiàn)，推薦閱讀Core et al., 2008.
而對(duì)于RNA-seq來(lái)說(shuō)，通過(guò)poly-A 等方法捕獲所得到的是成熟的mRNA，一方面，無(wú)法觀察到這些現(xiàn)象，另一方面所測(cè)到的成熟mRNA是細(xì)胞mRNA的累積量，而無(wú)法反應(yīng)某一時(shí)刻基因是否正在轉(zhuǎn)錄。
距離GRO-seq開(kāi)發(fā)出來(lái)，已經(jīng)過(guò)去了10余年。期間，測(cè)nascent transcription 的方法也在進(jìn)步和發(fā)展，例如和利用抗體富集的NET-seq（native elongating transcript sequencing）（Nojima et al., 2015）和利用biotin-labeled NTP 的PRO-seq（Mahat et al., 2016）可以實(shí)現(xiàn)單堿基分辨率的新生RNA鑒定; 可以精確鑒定TSS處RNA polymerase pausing site的PRO-Cap，ChRO-seq（Chu et al., 2018）。關(guān)于這些方法的詳細(xì)介紹和比較先留個(gè)坑，以后再填。今天我們主要想介紹一下GRO-seq的分析流程。
GRO-seq分析流程
同一般的高通量測(cè)序分析，GRO-seq的分析也主要分三步，首先進(jìn)行pre-processs reads，即去掉adaptor等序列；隨后進(jìn)行序列比對(duì)，最后分析比對(duì)結(jié)果。我們以mouse activated primary B cell 的GRO-seq (GEO: GSE116321,GSM3227964)為例進(jìn)行演示。一般來(lái)講，GRO-seq的建庫(kù)方法是單端測(cè)序的，讀長(zhǎng)較短為75bp。
預(yù)處理和序列比對(duì)
首先下載數(shù)據(jù)。使用ascp高速下載：
ascp -i ~/.aspera/connect/etc/asperaweb_id_dsa.openssh -Tr -Q \
anonftp@ftp-trace.ncbi.nlm.nih.gov:/sra/sra-instant/reads/ByRun/sra/SRR/SRR744/SRR7440347/SRR7440347.sra .
fastq-dump --split-3 --gzip SRR7440347.sra
然后進(jìn)行cutadapt trim adaptor
cutadapt -a GATCGGAAGAGCACACGTCT -a ACACTCTTTCCCTACACGACGCTCTTCCGATC -o SRR7440347.trimmed.fastq SRR7440347.fastq
使用bwa比對(duì)reads
filename=nsRNA_preB_GSM3227964
bwa aln -t 20 /home/hulab/chen/databases/bwa_indexes/mm10 SRR7440347.trimmed.fastq > ./SRR7440347.sai
bwa samse -n 1 /home/workdir/chen/databases/bwa_indexes/mm10 SRR7440347.sai SRR7440347.trimmed.fastq | samtools view -bS -q 30 - | samtools sort - > ${filename}.bam
samtools index ${filename}.bam
去除rRNA reads
這時(shí)，我們所得到的reads就是新生RNA鏈產(chǎn)生的了。但是注意哦，這里的reads不僅僅包括來(lái)自一般基因的RNA polymerase II 轉(zhuǎn)錄的，還包括了RNA polymerase I 轉(zhuǎn)錄的核糖體RNA，核糖體RNA轉(zhuǎn)錄的一般不是我們想要看的基因。需要對(duì)核糖體RNA進(jìn)行去除。
小鼠的核糖體rRNA可以從UCSC上獲得。
由于實(shí)驗(yàn)中已經(jīng)進(jìn)行了去除rRNA的步驟，因此去除rRNA 在總量上影響不大。
bed_rRNA=/home/workdir/chen/databases/mouse/genes_annotation/mm10_rRNA.bed
bedtools intersect -v -a ${filename}".bam" -b $bed_rRNA > ${filename}"_rm_rRNA.bam"
samtools view ${filename}"_rm_rRNA.bam" > ${filename}"_rm_rRNA.sam"
下游分析
下游的分析包括幾類。
針對(duì)于基因來(lái)講：新轉(zhuǎn)錄本的鑒定、基因的定量分析、promoter pausing index 的計(jì)算
和enhancer的分析，這里不展開(kāi)了.
使用homer做這些分析非常方便：
homer
首先先建立一個(gè)directory，利用makeUCSCfile 將bam轉(zhuǎn)換成bedGraph 格式，可以直接放在IGV里面可視化。
dir=$filename
###
makeTagDirectory $dir ${filename}"_rm_rRNA.sam" \
-genome mm10 -format sam -checkGC
# Estimated average read density = 0.004548 per bp
makeUCSCfile $dir/ -strand both -norm 1e6 -o ${filename}.bdg
然后對(duì)gene body 和TSS附近的reads 分別進(jìn)行定量分析。
#-------------- quantification of RNA and identify active genes ---------#
# the total reads are only 12 millions, so norm to 1e6
# the strand of gene is opposite to the reference.
analyzeRepeats.pl rna mm10 -count genes -condenseGenes -norm 1e6 -d $dir > ${filename}_genes_norm.txt
analyzeRepeats.pl rna mm10 -count pausing -condenseGenes -norm 1e6 -d $dir > ${filename}_pausing_norm.txt
結(jié)果分析
gene body 上的定量。 最后一列是normalize的 gene的RPM。
promoter 附近的定量。
這里定義的promoter 指 TSS 附近（-50，200）， 而gene body 所取的region是 （200，5000）。所謂Pausing ratio = TSS處的reads數(shù)/ gene body reads, 見(jiàn)倒數(shù)三列。這樣就知道了每個(gè)基因的附近 RNA polymerase的pausing 情況了。由于是比值，如果TSS read 和 gene body reads 都是0， 那么這個(gè)比值是1 ，這樣的情況應(yīng)該不被考慮。注意去除。
參考文獻(xiàn)
Core, L.J., Waterfall, J.J., and Lis, J.T. (2008). Nascent RNA Sequencing Reveals Widespread Pausing and Divergent Initiation at Human Promoters. Science 322, 1845–1848.
Nojima, T., Gomes, T., Grosso, A.R.F., Kimura, H., Dye, M.J., Dhir, S., Carmo-Fonseca, M., and Proudfoot, N.J. (2015). Mammalian NET-Seq Reveals Genome-wide Nascent Transcription Coupled to RNA Processing. Cell 161, 526–540.
Mahat, D.B., Kwak, H., Booth, G.T., Jonkers, I.H., Danko, C.G., Patel, R.K., Waters, C.T., Munson, K., Core, L.J., and Lis, J.T. (2016). Base-pair-resolution genome-wide mapping of active RNA polymerases using precision nuclear run-on (PRO-seq). Nature Protocols 11, 1455–1476.
Chu, T., Rice, E.J., Booth, G.T., Salamanca, H.H., Wang, Z., Core, L.J., Longo, S.L., Corona, R.J., Chin, L.S., Lis, J.T., et al. (2018). Chromatin run-on and sequencing maps the transcriptional regulatory landscape of glioblastoma multiforme. Nat Genet 50, 1553–1564.
pipeline from Danko Lab
https://github.com/Danko-Lab/tutorials/blob/master/PRO-seq.md
小tip：
GRO-seq使用BrU，研究的是以DNA為模板的轉(zhuǎn)錄情況。而加polyA 是在mRNA上直接加的，并且不會(huì)利用BrU。因此GRO-seq分析不用去除polydA 尾巴。
GRO-seq為什么使用短的讀長(zhǎng)：為了研究RNA polymerase 的移動(dòng)，降低了NTP的濃度，使得polymerase走的比較慢。最長(zhǎng)走80-100bp。
鏈特異性：是要在單鏈的時(shí)候就加上adaptor才能知道。如果是雙鏈加上adaptor則無(wú)法知道。
