RNA_seq pipline
RNA_seq pipline
PeRl
2018年3月7日
首先说明一下我做RNA-seq处理流程的文件树格式:
- RNA-seq/
- data/
- GRCh38.gtf
- chroms/
- hg38/
- samples/
- SraAccList.txt
- sra/
- fasta/
- fastqc/
- cufflinks_result/
- tophat_result/
- HTSeq_result/
- tools/
- Trimmomatic-0.36/
- data/
1. 下载参考基因组序列信息及注释文件GTF
参考基因组用于reads的定位和识别,需要用到的文件格式是fasta以.fa为后缀, 另外的各类文件存储格式可以参看file format.
chromFa.tar.gz是组装后的序列,每条染色体一个文件(我们要下载的文件),用axel下载数据
#例如我们要下载hg38,其中-n 20表示线程数
axel -n 20 http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chromFa.tar.gz解压后,我们需要把染色体文件整合成一个整体文件,然后用bowtie2建立索引文件
#! /usr/bin/env bash
for i in `ls /home/wang/Downloads/RNA-seq/data/chroms `
do
`cat /home/wang/Downloads/RNA-seq/data/chroms/$i >> /home/wang/Downloads/RNA-seq/data/hg38/hg38.fa`
done建立索引文件(有点耗时,约2小时左右):
bowtie2-build /home/test/RNA-seq/data/hg38/hg38.fa hg38从UCSC下载注释文件:
2..sra数据下载及格式转换
.sra数据的处理需要用到sra-toolkit
sudo apt-get install sra-toolkit然后用Bash脚本批量下载数据
#! /usr/bin/env bash
for i in `cat /home/wang/Downloads/RNA-seq/data/samples/SraAccList.txt`
do
echo $i
`axel -n 20 ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR212/$i/$i.sra > /home/wang/Downloads/RNA-seq/data/samples/sra/$i.sra`
# `mkdir /home/wang/Downloads/RNA-seq/data/samples/fasta/$i`
# `fastq-dump -O /home/wang/Downloads/RNA-seq/data/samples/fasta/$i --split-3 /home/wang/Downloads/RNA-seq/data/samples/sra/$i.sra`
done得到样本的.sra数据后,我们仍需要转换为.fa数据:
--split-3: 拆分文件,如果得到的.sra文件是单末端,那么这个参数就会被忽略;如果原文件中有两个文件,那么它就会把成对的文件按*_1.fastq,*_2.fastq这样分开。
fastq-dump [-O <outdir path>] –split-3 <input sra path>
fastq-dump -O /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441 --split-3 /home/wang/Downloads/RNA-seq/data/samples/sra/ERR2124441.sra3. Fasta质量控制
FastQC最为重要的部分是：Per base sequence quality（碱基质量）,Per base sequnence content（碱基含量分布）,Adapter Content（接头含量）。
sudo apt-get install fastqc
fastqc -q /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/*.fastq -o /home/wang/Downloads/RNA-seq/data/samples/fastqc- 总体的质量表现
上图中Summary的部分就是整个报告的目录,整个报告分成若干个部分。合格会有个绿色的对勾,警告是个“!”,不合格是个红色的叉子。
- 碱基质量
一般认为从第二个碱基开始,平均每个碱基的测序质量在四分位线在30分以上,则认为测序质量非常好,也就是我们通常称的Q30。
- 每条序列的测序质量统计
统计序列测序质量,如果均值在高分段就说明质量可以。
- 碱基含量统计
统计每种碱基含量,一般前面几bp序列可能包含引物而导致质量不高,可以cut掉；后面部分应该在25%左右,像下面的例子就是GC含量不合格。
- GC含量分布
蓝色的线是程序根据经验分布给出的理论值,红色是真实值,两个应该比较接近才比较好
- adapter
如果有adapter序列没有去除干净的情况,在后续分析的时候需要先使用Trimmomatic软件。
Trimmomatic需要区分PE(Paired End)和SE(Single End).
java -jar <path to trimmomatic.jar> PE/SE [-threads <threads] [-phred33 | -phred64] [-trimlog <logFile>] <input 1> <input 2> <paired output 1> <unpaired output 1> <paired output 2> <unpaired output 2> <step 1> …
过滤数据的步骤与命令行中过滤参数的顺序有关,参数之间用空格分隔,参数和参数数值之间用:分隔,通常的过滤步骤如下：
| 参数 | 意义 | options |
|---|---|---|
| ILLUMINACLIP | 过滤reads中的Illumina测序接头和引物序列,并决定是否去除反向互补的R1/R2中的R2 | <需要的adapter地址>:<最大错配数>:<跨adapter序列的匹配精确度>:<两个adapter之间序列的匹配精确度> |
| SLIDINGWINDOW | 从 reads 的 5’ 端开始,进行滑窗质量过滤,切掉碱基质量平均值低于阈值的滑窗 | <滑动窗大小>:<质量平均值> |
| MAXINFO | 一个自动调整的过滤选项,在保证 reads 长度的情况下尽量降低测序错误率,最大化 reads 的使用价值 | NULL |
| LEADING | 从 reads 的开头切除质量值低于阈值的碱基 | <最小质量> |
| TRAILING | 从 reads 的末尾开始切除质量值低于阈值的碱基 | <最小质量> |
| CROP | 从 reads 的末尾切掉部分碱基使得 reads 达到指定长度 | <长度> |
| HEADCROP | 从 reads 的开头切掉指定数量的碱基 | <碱基数量> |
| MINLEN | 如果经过剪切后 reads 的长度低于阈值则丢弃这条 reads | <最短长度> |
| AVGQUAL | 如果 reads 的平均碱基质量值低于阈值则丢弃这条 reads | NULL |
| TOPHRED33 | 将 reads 的碱基质量值体系转为 phred-33 | NULL |
| TOPHRED64 | 将 reads 的碱基质量值体系转为 phred-64 | NULL |
java -jar /home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/trimmomatic-0.36.jar PE -threads 12 -phred33 -trimlog /home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/logfile.txt /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1.fastq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2.fastq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1_paired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_1_unpaired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2_paired.fq /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/ERR2124441_2_unpaired.fq ILLUMINACLIP:/home/wang/Downloads/RNA-seq/tools/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:504. Tophat mapping
处理完reads后，我们需要用tophat把reads匹配到参考基因组上。
tophat -o <output path> -G <gtf path> <索引文件path> <fasta文件path>
tophat -o /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441 -G /home/wang/Downloads/RNA-seq/data/GRCh38.gtf /home/wang/Downloads/RNA-seq/data/hg38/hg38 /home/wang/Downloads/RNA-seq/data/samples/fasta/ERR2124441/*_paired.fq最后得到的结果如下，一般情况下，接下去需要用到的就是accepted_hits.bam
接着我们需要检查一下得到的accepted_hits.bam是否已经排序:
sudo apt-get install samtools
samtools view -h /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam | head如果header中显示SO:coordinate,那就说明已经排序了。
5. RNA-Seq 质控
基于原始测序数据的质控（例如FastQC）不足以保证RNA-seq 数据的可用性。
我们还需要使用RSeQC程序包(基于python)全面地评估RNA-seq 结果质量，例如序列质量、GC 偏倚、PCR 偏倚、核苷酸组成偏倚、序列深度、链特异性、覆盖均一性，和基因组结构上的片段分布等，以确保后续分析的可靠性。
#安装RSeQc
sudo apt-get install python-pip
pip install RSeQc
#如果出现compr_lzo.c:31:23: fatal error: lzo/lzo1x.h: No such file or directory,需要安装 liblzo2-dev
sudo apt-get install liblzo2-dev简单的用bam_stat.py看下统计信息
[bam_stat.py的路径] -i [.bam的路径] > [输出路径]
#找到python安装包的文件夹地址
.local/bin/bam_stat.py -i /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam > /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/output_bam_stat.txt
6. cufflinks得到fpkm值
因为前面用的GRch38.gtf用的是Refseq ID,最后还需要ID转换，这个就不赘述了。
#sudo apt-get install cufflinks
cufflinks -p 4 -G /home/wang/Downloads/RNA-seq/data/GRCh38.gtf -o /home/wang/Downloads/RNA-seq/data/samples/cufflinks_result/ERR2124441 /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam7. HTSeq得到count数及DEseq2处理
处理count数用到的是python工具HTSeq
pip install HTSeqhtseq-count使用时需用samtools对accepted_hist.bam文件以read name排序（-n）才能使用
samtools sort -n /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits.bam -o /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits_sort.bam 排序后用htseq-count计数,输入文件为排序后的tophat比对结果accepted_hits_sort.bam
根据不同需要可以修改参数，参看官方说明,其中-i gene_id描述的是使用的:
/home/wang/.local/bin/htseq-count -f bam -s no -a 10 -t exon -i gene_id -m union /home/wang/Downloads/RNA-seq/data/samples/tophat_result/ERR2124441/accepted_hits_sort.bam /home/wang/Downloads/RNA-seq/data/samples/cufflinks_result/ERR2124441/transcripts.gtf > /home/wang/Downloads/RNA-seq/data/samples/HTSeq_result/ERR2124441/SRR1206035.htseq_count
接下去的处理就是用DESeq2处理count数，以前我有做过总结，这方面就不再加以说明了，就是酱紫~\(≧▽≦)/~