Ribo-seq 数据上游分析

点击上方关注老俊俊

1引言

看过上面推文的小伙伴可能不知道上面的 bam 文件从何而来,虽然我在文中有描述过,也有小伙伴问我,这次介绍一下上游数据的处理。

2步骤

大概步骤如下:

1. fastqc 质控。
2. 去接头(optional,根据自己接头序列选择)。
3. 去除 rRNA/tRNA 序列。
4. 比对到(基因组/转录本序列)。
5. 定量及下游分析。

3qc 质控

使用 fastqc 和 multiqc 对原始测序数据进行 qc 检测和整合,这里给个示例代码:

# qc
$ fastqc -t 10 -q -o 0.QC-data/raw-qc/ 1.raw-data/*.gz

# combine
$ multiqc 0.QC-data/raw-qc/* -o 0.QC-data/raw-qc/

4去接头

这里根据自己的接头或者文章里的接头序列进行去除,软件也很多,cutadapt, trim_galore 等:

# trim adapters
$ trim_galore -j 5 -q 20 \
              --no_report_file \
              --length 20 --max_length 35 \
              -o 2.trim-data/ \
              1.raw-data/test.fq.gz

5去除 rRNA/tRNA

因为核糖体保护的片段很短,大约在 29nt 左右,在提取过程中会包含大量的其它序列的污染,而 rRNA 和 tRNA 序列在细胞里是站绝大部分的(rRNA),为了避免这样的污染, 防止测序过程中扩增产生大量这些无用的数据, 一般会在建库过程中使用 rRNA 试剂盒来进行去除,如果没有做的话,或者在测序数据拿到以后将这些比对到 rRNA/tRNA 序列上进行去除,拿到有用的数据。

NCBI 下载对应物种的 rRNA/tRNA 序列

---

同样的方法把 tRNA 序列也下载下来,重命名一下即可:

建立索引

# rRNA索引
$ mkdir -p index/rRNA-index
$ bowtie2-build --threads 20 human_rRNA.fasta -p index/rRNA-index/human_rRNA

# tRNA索引
$ mkdir -p index/tRNA-index
$ bowtie2-build --threads 20 human_tRNA.fasta -p index/tRNA-index/human_tRNA

比对保留未比对的序列

# rm rRNA
$ mkdir 3.rmrRNA-data
$ bowtie2 -p 22 -x index/rRNA-index/human_rRNA \
          --un-gz 3.rmrRNA-data/test.rmrRNA.fq.gz \
          -U 2.trim-data/test-trim.fq.gz \
          -S 3.rmrRNA-data/test.sam \
          && rm 3.rmrRNA-data/test.sam

# rm tRNA
$ mkdir 4.rmtRNA-data
$ bowtie2 -p 22 -x index/tRNA-index/human_tRNA \
          --un-gz 4.rmtRNA-data/test.rmtRNA.fq.gz \
          -U 3.rmrRNA-data/test.rmrRNA.fq.gz \
          -S 4.rmtRNA-data/test.sam \
          && rm 4.rmtRNA-data/test.sam

根据经验,一般比对到 rRNA 上的比对率在 70-90 左右,比对到 tRNA 上的比对率在 30 左右,当然也和实验有关,只是一个大概。

6比对到基因组/转录组

如果要得到开头提到的教程里的类似的数据格式,我们得比对到转录组序列上面,此外转录组序列也是需要筛选过的。

筛选的代码我已经写成软件叫做 GetTransTool,见推文 python package: GetTransTool。

安装

$ pip install GetTransTool

提取代表转录本序列

这里使用 GetCDSLongestFromGTF 函数从GTF文件及基因组文件提取出 CDS 区最长 的代表转录本序列:

$ GetCDSLongestFromGTF --database ensembl \
                       --gtffile Homo_sapiens.GRCh38.101.gtf.gz \
                       --genome Homo_sapiens.GRCh38.dna.primary_assembly.fa \
                       --outfile longest_cds_trans.fa

Your job is running, please wait...
Your job is done!
Running with 152.38 seconds!

提取到的转录本序列如下:

>OR4F5_ENSG00000186092_ENST00000641515_61_1038_2618
CCCAGATCTCTTCAGTTTTTATGCCTCATTCTGTGAAAATTGCTGTAGTCTCTTCCAGTT
ATGAAGAAGGTAACTGCAGAGGCTATTTCCTGGAATGAATCAACGAGTGAAACGAATAAC
TCTATGGTGACTGAATTCATTTTTCTGGGTCTCTCTGATTCTCAGGAACTCCAGACCTTC
CTATTTATGTTGTTTTTTGTATTCTATGGAGGAATCGTGTTTGGAAACCTTCTTATTGTC
...

格式为 基因名-基因id-转录本id-CDS起始位置-CDS终止位置-转录本长度 。

建索引

# 转录组索引
$ mkdir -p index/trans-index
$ bowtie-build index/longest_cds_trans.fa -p index/trans-index/trans

比对到转录组

前面所需东西准备好了即可比对:

注意 bowtie 所需输入文件必须是 非压缩 的 fastq 文件!

$ zless 4.rmtRNA-data/test.rmtRNA.fq.gz | \
  bowtie index/trans-index/trans \
         -q -p 20  -m 1 -v 2 \
         --best --strata - -S \
         |samtools sort -@ 20 -o 5.map-data/test.sorted.bam

# index bam
$ ls 5.map-data/*.bam |xargs -i samtools index -@ 20 {}

查看 bam 为文件:

$ samtools view test.sorted.bam |less -S
SRR548.8675742  0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       738     255     30M     *       0  >
SRR548.165589   0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       799     255     33M     *       0  >
SRR548.3116679  0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       804     255     21M     *       0  >
SRR548.11512505 0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       804     255     28M     *       0  >
SRR548.13404411 0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       807     255     28M     *       0  >
SRR548.12704826 0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       829     255     33M     *       0  >
SRR548.3045804  0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       842     255     31M     *       0  >
SRR548.9179985  0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       847     255     26M     *       0  >
SRR548.763651   0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       899     255     29M     *       0  >
SRR548.10238203 0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       961     255     33M     *       0  >
SRR548.5980239  0       SDF4_ENSG00000078808.19_ENST00000360001.12_SDF4-202_285_1352_1933       969     255     31M     *       0
...

7结尾

这里的 bam 文件即可用开头的推文里的代码进行可视化了,此外需要注意位置计算问题(frame 分配),大家可以自己去检查尝试。这里就不废话了。

欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 哦,数据代码已上传至QQ群,欢迎加入下载。

群二维码：

老俊俊微信：

知识星球：

所以今天你学习了吗？

今天的分享就到这里了，敬请期待下一篇！

最后欢迎大家分享转发，您的点赞是对我的鼓励和肯定！

如果觉得对您帮助很大，赏杯快乐水喝喝吧！

 往期回顾 

◀...