RNA-seq分析(1)之数据质控

本文介绍RNA-seq上游分析之数据质控：

fastqc质控
- 去除接头序列
- reads trim
去除rRNA序列
- rRNA index
- 反向mapping

1.fastq文件质控

1	fastq -f fastq -t 10 -o ./ ./myfastq_r1.fq.gz ./myfastq_r2.fq.gz

之后会产生质控报告和zip文件，查看报告，可以看到base quality、GC content等信息。质控报告每个版块的具体的意义参考(https://zhuanlan.zhihu.com/p/20731723 )。

2.去除接头序列

RNA-seq建库测序时，会在样品末端加上接头序列，在后续的步骤中需要去除。

Usage:
    cutadapt -a ADAPTER [options] [-o output.fastq] input.fastq

For paired-end reads:
    cutadapt -a ADAPT1 -A ADAPT2 [options] -o out1.fastq -p out2.fastq in1.fastq                            in2.fastq

去除接头序列，这里我用cutadapt，对于双端测序来说，需要去掉read1和read2两个序列，参数-a和-A分别是read1和read2序列的3’端序列。-e即error rate,-n为2即最多每条read最多去掉两个adapter，-m为20去掉adapter之后若reads长度小于20，则去掉这条reads，–pair-filter这个参数对于双端测序而言，read1和read2都有可能检测到接头。如果选择any，则只要两个中其中一个检测到接头，read1和read2均舍弃；如果选择both，则必须两个都检测到接头，read1和read2才舍弃，-j线程数，-q [5’CUTOFF,]3’CUTOFF参数，quality cutoff，去掉5’端、3’端碱基质量低于该值。双端测序需要提供两个值。

cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -A AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT
\ -e 0.1 -O 5 -n 2 -m 20 --pair-filter both -q 20,20 -j 10
\ -o 01.clean/CS-1.clean.r1.fq.gz -p 01.clean/CS-1.clean.r2.fq.gz
\ 00.rawdata/CS-1-LEG8830_L2_1.fq.gz 00.rawdata/CS-1-LEG8830_L2_2.fq.gz

3.去除测序刚开始时不稳定的reads

测序仪进行测序时，刚开始10bp左右的reads测序不稳定，可以进行适当修剪。

1	cutadapt -u 10 -o trimmed.fastq input_reads.fastq

4.去除rRNA序列

有时候质控得到的结果中GC含量为双峰，这可能是一些RNA过表达引起的GC偏好，需要去除一些rRNA序列或是其他我们并不关心的序列，但是也并不必过于担忧，RNA-seq中本身就有一些表达量高的序列，所以质控报告中只要一些如碱基质量等的关键因素满足要求，后续只是进行差异表达等的分析，不去除也是可以的。
首先，我们需要在NCBI数据库中下载rRNA序列的fasta文件，然后用Hisat2进行比对，保留没有比对上的序列即可,过程如下：

4.1 构建索引

1	hisat2-build -p 10 rRNA.fasta ./rRNA/rRNA

4.2 输出没有比对上的序列

–un-conc-gz参数输出没有比对上的序列。

for i in CS-1 CS-2 PS-1 PS-2
do
	hisat2 -x 03.rm/rRNA/rRNA \
	-1 01.clean/$i.r1.clean.fq.gz \
	-2 01.clean/$i.r2.clean.fq.gz \
	--un-conc-gz 03.rm/${i}_rmr_%.fq.gz -p 10 -S 03.rm/${i}.rRNA.sam
done

这一步进行完之后便会发现fastq文件小了许多，这便是去除成功了。然后便可以继续接下来的分析了。

参考资料：
1.https://zhuanlan.zhihu.com/p/20731723
2.https://mp.weixin.qq.com/s/jVTmp4VJQxB6Lwxh1QmNkA