今天开始正式进入RNA-seq分析啦~
软件安装
由于昨天就已经创建好了新的conda的环境，所以今天只要把用到的软件都装一下就好了。
把要安装的软件分一下类型：
数据格式转换
-
sra-tools: 这个软件用于把NCBI里下载的.sra格式的数据给转换成fastq。用到的工具是里面的fastq-dump
数据质控
数据质量的好坏是数据分析中的关键，所谓“garbage in, garbage out”，如果一开始的数据是有问题的，后续再怎么分析都是徒劳。
那数据会存在什么样的问题呢？
主要是二代测序增加进去adapters序列，和测序中的低质量数据。它们会影响后续的比对、定量准确性，更会影响下游分析，因此要先生成质量报告了解一下数据的质量，如果质量不好的话要把不符合要求的序列给过滤掉。
所以这一块分两步，1. 生成质量报告；2. 过滤不符合要求的序列。
-
fastqc&multiqc: 用于生成质量报告。fastqc能对每个样品生成单独的质量报告，而multiqc则是可以把多个样本的fastqc结果整合到一个文件中，这样比较方便全面的了解和对比多个数据之间的质量情况。 -
fastp: 由陈实富大佬写的新手友好型质控软件，只要设置输入文件和输出文件就好了，推荐给各位胖友。之前写过一篇关于这个软件的简书文章：用fastp对转录组数据做QC -
trim_galore&trimmomatic另外两个用于做质控的软件。为了防止质控软件本身的算法原因导致分析结果的误差，条件允许的情况下建议用不同的质控软件来处理数据进行对比，取结果好的进行下游分析。
比对（mapping）
有参转录组和无参转录组分析的本质区别就在于有没有这个“参”。参是什么呢？是reference，中文世界一般叫“参考序列”。而我研究生期间做的植物基因组，就是在做这个“参”。
比对的软件有很多种，有基于基因组比对（star、hisat2）、基于转录组比对(bowtie、bwa)、不基于比对(salmon)
我目前用bwa比较多，别的只是知道，没怎么用过。
定量
-
subread用的是其下的featureCounts这个小程序。这一块我不太了解，这也是我参加这次转录组学习小组的原因~
用conda安装软件
# 激活专门用于RNA-seq分析的conda环境
conda activate rnaseq
# 安装软件
conda install fastp fastqc multiqc trim-galore trimmomatic hisat2 bowtie2 subread salmon
因为我目前人在Singapore，所以用conda下载软件还挺快的。有的小朋友习惯在conda install的时候加-y参数，这样就只要等着它装好就好了，不需要再手动去确定是否安装这些软件。但是也有个问题，它在整个安装过程中不会给任何的提示和输出，就只能看到solving environment 的状态栏在“爱的魔力转圈圈”，当安装很多软件或者网络不好的时候会特别没有安全感。
数据格式转换
我想想我一般是怎么写这个的：
#!/usr/bin/env bash
list="SRR1039509.sra
SRR1039512.sra
SRR1039513.sra
SRR1039516.sra
SRR1039517.sra
SRR1039520.sra
SRR1039521.sra"
for i in $list
do
fastq-dump --gzip --split-3 ${i} -O ../01raw
done
我的list是通过下面的命令生成的:
ls *.sra | xargs -l
而豆豆提供的脚本是这个画风：
raw=~/rnaseq/raw # 这个路径根据个人需求进行修改即可
cat $raw/srr.ids | while read i ;do
echo $i
# time fastq-dump --gzip --split-3 -A $i $raw/${i}.sra -O $raw;
done
确实看起来比我高明得多呢。学习了学习了。用time这个命令计算运行时间是一个很好的习惯。
先用井号注释掉了真正要运行的命令，这在写脚本的过程中是个好习惯，因为有时自己会搞错变量(比如这里的i) ，于是用echo $i可以先看看是不是真正要用的，先运行一下看看结果
把脚本写进一个叫做fastqdump.sh的文件，然后nohup运行。
(不要随便起个名字噢，不然过段时间就忘了这个当时写了是做什么用的了呢，后患无穷，整理起来特别头疼。)
我用nohup主要是它可以保留软件运行过程中的日志和报错，当然不用nohup也可以做到这点：
bash fastqdump.sh 1>fastqdump.log 2>&1
# 或者
bash fastqdump.sh 1>fastqdump.log 2>fastqdump.err
1>fastqdump.log 表示将结果的正确日志文件输出到fastqdump.log中，2>&1 这个2表示错误日志，将它也合并到1的正确日志中。
当然，如果不想要错误信息，可以把它丢到linux"黑洞"中去
2>/dev/null
下载参考序列
根据文章介绍，所用的reference是hg19的版本。因此就下载这个版本的genome和gtf数据备用。
wget -c http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.fa.gz
wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
太晚了。明天继续写。
Hans Chen
2019年6月9日于 Westwood Hostel, Jurong West Ave. 5, Singapore.
网友评论