数据的准备
RNA-seq主要是用来检测不同的时空或者不同的状态下的基因表达差异。常用流程有topHat+Stringtie+Ballgown。这里将主要使用Hisat2+featureCount+DEseq2的流程(同时参照其他流程把其他软件选择加入),可能更(gèng)新。
准备原始数据
本次RNA-seq流程来自Github的twbattaglia,基本照着以上的流程进行。
下载注释文件以及基因组。目前最新的小鼠(家鼠)参考基因组是GRCm39(mm39),但是可能还是GRCm38(mm10)用得比较多。本着向前发展的原则,我们用mm39。
| wget -P anntation/ http://hgdownload.soe.ucsc.edu/goldenPath/mm39/bigZips/genes/mm39.ncbiRefSeq.gtf.gz
wget -P genome/ http://hgdownload.soe.ucsc.edu/goldenPath/mm39/bigZips/mm39.fa.gz
gunzip anntation/mm39.ncbiRefSeq.gtf.gz
gunzip genome/mm39.fa.gz
|
然后是数据,数据来源是来自“在低葡萄糖和高葡萄糖条件下培养的小鼠胰岛的转录组差异分析”。文章PMID25051960
| wget -P rawdata/ ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR137/001/SRR1374921/SRR1374921.fastq.gz
wget -P rawdata/ ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR137/002/SRR1374922/SRR1374922.fastq.gz
wget -P rawdata/ ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR137/003/SRR1374923/SRR1374923.fastq.gz
wget -P rawdata/ ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR137/004/SRR1374924/SRR1374924.fastq.gz
|
建立索引
使用hisat2对小鼠基因组mm39.fa建立索引。
| hisat2-build genome/mm39.fa mm39
|
使用STAR对小鼠基因组mm39.fa建立索引,建立的同时把gtf加入了,比对后应该可以直接输出counts统计。
| STAR \
--runMode genomeGenerate \
--genomeDir genome/star_index \
--genomeFastaFiles genome/mm39.fa \
--sjdbGTFfile anntation/mm39.ncbiRefSeq.gtf \
--runThreadN 16
|
我对STAR没啥好感,因为它的--help又臭又长。另外跑个建立索引就把我的小机器的16G内存吃完了,但是CPU占用率又上不去。
质控
对RNA-seq进行质控的操作与DNA-seq的质控基本一样,使用FastQC等软件来查看数据质量,再用trimmomatic或trim_galore等软件来进行数据过滤。
单端数据就这样跑好了,双端也是类似的。
fastqc基本用法
| fastqc rawdata/SRR1374921.fastq.gz \
-o QC/ --noextract
|
上述命令获得质量报告后,使用trim galore进行过滤。注意,trim galore必须在cutadapt已在环境变量中才可以使用。
| trim_galore \
--quality 20 \
--fastqc \
--length 25 -j 10 \
--output_dir cleandata \
rawdata/SRR1374921.fastq.gz
|
对所有样本均进行以上操作。
把这些样本ID都扔进一个samplelist里就可以批量运行了。
| cat samplelist.txt | while read line;
do
fastqc rawdata/${line}.fastq.gz -o QC/ --noextract;
trim_galore --quality 20 --fastqc --length 25 -j 10 \
--output_dir cleandata rawdata/${line}.fastq.gz;
done
|
去除rRNA(可选)
接下来使用SortMeRNA去除rRNA,这一步其他的RNA-seq分析流程比较少见。
从sortmerna的原始官网中,是可以下载到教程中的2.1版本的,同时下载下来里面也有silva的16s等数据库存在。但在更新的github 4.2版本中,提供的数据库已经不是2.1版本中提供的了,所以这里接下来将不按原教程的操作,而是使用4.2版本进行。
使用的数据库为
| wget https://pd.zwc365.com/seturl/https://github.com/biocore/sortmerna/releases/download/v4.2.0/database.tar.gz
tar -zxvf database.tar.gz
|
使用里面的smr_v4.3_default_db.fasta。
| sortmerna --ref database/smr_v4.3_default_db.fasta \
--reads cleandata/SRR1374921_trimmed.fq.gz \
--aligned rRNA/SRR1374921_aligned \
--other rRNA/SRR1374921_filtered \
--fastx -threads 10 -v \
--workdir temp/
gzip rRNA/SRR1374921_aligned.fq
gzip rRNA/SRR1374921_filtered.fq
rm -rf temp/
|
这个软件使用会自动创建Key-value DB和索引文件夹,后面加上--workdir参数来指定一个临时文件夹来存放。另外这软件虽然能输入gz压缩文件,但是输出只能是fastq。得到的filtered.fq将在后续进行使用。