实践RNA-seq

安装软件

安装miniconda

进入清华大学镜像站内下载一个合适的版本

wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh
bash Miniconda3-latest-Linux-x86_64.sh
echo $PATH #查看环境变量

配置镜像

conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/pkgs/free
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/conda-forge
conda config --add channels https://mirrors.tuna.tsinghua.edu.cn/anaconda/cloud/bioconda
conda config --set show_channel_urls yes

创建名为rna的软件安装环境，同时安装python2版本

conda create -n rna python=2

激活/进入conda的rna环境

source activate rna

安装 sra-tools软件

conda install -y sra-tools # 接下来要用它下载sra数据，这是sra数据所在网页上要求安装的。

数据下载

根据文献找到GSE号，这里是GSE52778

airway文献.png

把文章中提到的GSE编号（GSE52778）提交到GEO数据库。https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52778
在页面中找到SRA后面的链接，即点击SRP033351字样，在弹出的页面中点击字样Send results to Run selector，在弹出的页面点击Accession List即可得到下载得到SRR_Acc_List.txt，下载这个list里的文件。文件默认下载到~/ncbi/public/sra

cd ~
mkdir -p rna.GSE52778/sra
cd rna.GSE52778/sra
source activate rna
cat SRR_Acc_List.txt | while read id; do (prefetch ${id});done &

数据下载的提速设置使用 aspera

以下内容来此 AmyCui_48c4

wget -c https://download.asperasoft.com/download/sw/connect/3.8.1/ibm-aspera-connect-3.8.1.161274-linux-g2.12-64.tar.gz
tar zxvf aspera-connect-3.7.4.147727-linux-64.tar.gz
bash aspera-connect-3.7.4.147727-linux-64.sh
cd # 去根目录
ls -a # 如果看到.aspera文件夹，代表安装成功
 
# 永久添加环境变量
echo 'export PATH=~/.aspera/connect/bin:$PATH' >> ~/.bashrc
# 最好写绝对路径：echo 'export PATH=/home/qmcui/.aspera/connect/bin:$PATH' >> ~/.bashrc
source ~/.bashrc
# 查看帮助文档
ascp --help
# 再用prefetch下载数据的时候，就不会调用https  而是使用fasp

数据转换（sra转fastq）

# 为了目录整洁，我把sra文件放到 ~/rna.GSE52778/sra
nohup ls /four/zy/rna.GSE52778/sra/*sra | while read id; do (fastq-dump --gzip --split-3 -O /four/zy/rna.GSE52778/sra.fq ${id}); done 1>sra2fq 2>&1 &

质控

data statistics

# 软件没有的话，就安装一下
ls /four/zy/rna.GSE52778/sra.fq/*gz |while read id; do (fastqc -t 2 -o /four/zy/rna.GSE52778/fq.qc ${id}); done

整合结果

# 软件没有的话，就安装一下
multiqc ~/rna.GSE52778/fq.qc

filter data

前一个步骤的结果为这一步骤的参数设置提供了依据。先建立一个config文件

ls ~/rna.GSE52778/sra.fq/*_1.fastq.gz >1
ls ~/rna.GSE52778/sra.fq/*_2.fastq.gz >2
paste 1 2  > config

考虑到希望为每个样品都生成一个日志信息，那么就先要为config文件内容再增加一列，这一列是对应的样品名称，这样共3列，由于是3列内容，我将之命名为config3。另外在rna.GSE52778之下新建一个叫做clean的路径，这一步的处理结果都放这儿。

paste config ~/rna.GSE52778/sra/SRR_Acc_List.txt >config3
cat config3 #看一下是不是想要的内容。至此，clean目录下有一些临时文件，比如1，2，config，删除不要了。

再建立一个bash脚本文件，命名为qc.sh，vim写入内容如下

cat $1|while read id
do
    arr=(${id})
    fq1=${arr[0]}
    fq2=${arr[1]}
    sample=${arr[2]}
    echo "sample ID is " $sample
    echo "################################################" 
    trim_galore -q 25 --phred33 --length 36  --stringency 3 --paired -o ~/rna.GSE52778/clean \
    $fq1 $fq2  1>$sample.trim.log 2>&1 # 不写日志信息也可以，该软件也自动生成相关信息。
    echo "OK"
done

这个脚本的建立也可以cat >qc.sh 然后贴入提前写好的内容，Ctrl C退出即可。最后cat qc.sh确认一下内容。运行名为qc.sh脚本

nohup bash qc.sh config3 &

比对

将reads比回参考基因组，得到比对的结果为sam文本，用于后续分析。sam文件紧接着转为bam文件。软件使用hisat2，参考基因组为 hg38.fa。先用hisat2为 hg38.fa构建索引。自己构建太慢，可以下载到。hisat2官网下载

wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/grch38.tar.gz # 下载之后解压

hisat2比对到索引文件，然后samtools直接对生成的bam文件排序，以利于后续软件分析。

复制一份SRR_Acc_List.txt到clean路径里，在此路径下执行：

nohup cat SRR_Acc_List.txt |while read id;do
hisat2 -p 5 -x ~/index/grch38/genome -1 \
${id}_1_val_1.fq.gz -2 ${id}_2_val_2.fq.gz | \
samtools sort -@ 5 -o ~/rna.GSE52778/sort.bam/${id}.sort.bam -
done &

差异分析

使用featureCounts，该软件包含在subread软件内。

# 在GSE52778_RNA下面建立featureCounts路径
# 确定已经下载好了基因组的gtf文件。这是在家目录下的gtf路径里
nohup featureCounts -T 5 -p -t exon -g gene_id \
-a ~/gtf/gencode.v29.annotation.gtf \
-o ~/rna.GSE52778/featureCounts/all.counts.txt SRR1039*.sort.bam &