Ribo-seq 数据上游分析

测试开头

点击上方关注老俊俊

1引言

Ribo-seq 数据质控研究-上

Ribo-seq 数据质控研究-中

Ribo-seq 数据质控研究-下

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

2步骤

大概步骤如下:

fastqc 质控。

去接头(optional,根据自己接头序列选择)。

去除 rRNA/tRNA 序列。

比对到(基因组/转录本序列)。

定量及下游分析。

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
...