RNA测序(RNA-seq)是在next-generation sequencing发明以来很快被人们肯定的一种研究生物表达组学的分析手段。它的出现深远地影响着2005年以后的基因表达与调控的研究,它提供了一个全面筛选目的基因全新手段。其基本原理如下图所示:
RNA-seq protocol
在本文中,我们不重点介绍RNA-seq相关原理,而仅关注其上游处理和分析步骤,目前对于RNA-seq数据的处理流程很多 但是都大致包括:质控、去除接头和低质量序列、maping、count的计算等等。
本文需要用到的软件有:
SRA Toolkit https://www.ncbi.nlm.nih.gov/sra/docs/toolkitsoft/
fastqc https://www.bioinformatics.babraham.ac.uk/projects/fastqc/
multiqc https://multiqc.info
HISAT2 http://ccb.jhu.edu/software/hisat2
samtools http://samtools.sourceforge.net/
StringTie http://ccb.jhu.edu/software/stringtie
首先我们使用 SRA Toolkit工具包中的prefetch工具下载GEO中的数据,并且使用fastq-dump将其释放为fastq文件
mkdir rnaseq prefetch SRR8382645 fastq-dump --gzip --split-3 SRR8382645 #但是由于fastq-dump只能利用单线程进行计算,速度实在不太理想,因此我们考虑用更为快速的fasterq-dump替代它 fasterq-dump --split-3 SRR8382645 #如果你是用的后者的话,稍后可以压缩产生的文件 gzip SRR8382645.sra_1.fastq
接着,我们使用fastqc对数据进行质控,并使用MultiQC整合fastqc的结果
mkdir fastqc_res fastqc -o fastqc_res SRR*.fastq.gz multiqc fastqc_res
使用TrimGalore去除adaptor
trim_galore --paired -o output/ SRR8382645.sra_1.fastq.gz SRR8382645.sra_1.fastq.gz
Mapping:
接下载使用hisat2进行maping
命令格式参考如下
hisat2 -p 8 --dta -x indexes/hg19/genome -1 \ /removed_trimed/SRR8382645.sra_1_trimed.fastq.gz \ -2 /removed_trimed/SRR8382645.sra_2_trimed.fastq.gz -S output/SRR8382645.sam
采用samtools进行排序和转换
samtools sort [email protected] 8 -o output/SRR8382645.bam output/SRR8382645.sam
注意该命令可能不适合 1.2.1 版本以下的samtools,在这些版本上可能需要先转换再进行排序,例如:
samtools view -bS output/SRR8382645.sam >output/SRR8382645_unsorted.bam
排序操作如下
samtools sort [email protected] 8 output/SRR8382645.bam SRR8382645
未完待续…