Skip to content

数据的准备

RNA-seq主要是用来检测不同的时空或者不同的状态下的基因表达差异。常用流程有topHat+Stringtie+Ballgown。这里将主要使用Hisat2+featureCount+DEseq2的流程（同时参照其他流程把其他软件选择加入），可能更（gèng）新。

准备原始数据

本次RNA-seq流程来自Github的twbattaglia，基本照着以上的流程进行。

下载注释文件以及基因组。目前最新的小鼠（家鼠）参考基因组是GRCm39（mm39），但是可能还是GRCm38（mm10）用得比较多。本着向前发展的原则，我们用mm39。 

1
2
3
4
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 

1
2
3
4
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建立索引。

1
hisat2-build genome/mm39.fa mm39

使用STAR对小鼠基因组mm39.fa建立索引，建立的同时把gtf加入了，比对后应该可以直接输出counts统计。

1
2
3
4
5
6
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基本用法 

1
2
fastqc rawdata/SRR1374921.fastq.gz \
    -o QC/ --noextract

上述命令获得质量报告后，使用trim galore进行过滤。注意，trim galore必须在cutadapt已在环境变量中才可以使用。

1
2
3
4
5
6
trim_galore \
    --quality 20 \
    --fastqc \
    --length 25 -j 10 \
    --output_dir cleandata \
    rawdata/SRR1374921.fastq.gz

对所有样本均进行以上操作。

把这些样本ID都扔进一个samplelist里就可以批量运行了。 

1
2
3
4
5
6
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版本进行。

使用的数据库为 

1
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。

1
2
3
4
5
6
7
8
9
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将在后续进行使用。