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.sra

3. 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> …

过滤数据的步骤与命令行中过滤参数的顺序有关,参数之间用空格分隔,参数和参数数值之间用:分隔,通常的过滤步骤如下：

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:50

4. 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.bam

7. HTSeq得到count数及DEseq2处理

处理count数用到的是python工具HTSeq

pip install HTSeq

htseq-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数，以前我有做过总结，这方面就不再加以说明了，就是酱紫~(≧▽≦)/~