安装软件
安装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 &
网友评论