RNA-seq上游分析流程

在本次pipeline中所使用的软件为：fastqc -- multiqc -- trim_galore -- rsem-STAR

一、conda安装

首先要进行conda的安装和环境配置

①下载conda

wget -c https://repo.continuum.io/miniconda/Miniconda3-latest-Linux-x86_64.sh
#更改权限
chmod 777 Miniconda3-latest-Linux-x86_64.sh

②安装和激活conda

bash Miniconda3-latest-Linux-x86_64.sh
export PATH="/home/sysbio/SA20008254/y/bin:$PATH" 
source ~/y/miniconda3/bin/activate 

③配置conda镜像源，由于目录中并无全局~/.condarc,因此要创建.condarc文件。

#创建默认镜像源
conda config --add  channels  defaults
#添加镜像源，用北外的（清华速度较慢）
conda config --add channels http://mirrors.bfsu.edu.cn/anaconda/cloud/bioconda/
conda config --add channels http://mirrors.bfsu.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels http://mirrors.bfsu.edu.cn/anaconda/pkgs/free/
conda config --add channels http://mirrors.bfsu.edu.cn/anaconda/pkgs/main/

④配置conda环境，这里跑RNA-seq，因此环境名字设置为RNA-seq，之后选择y。

conda create -n RNA-seq python=3.8
#激活环境
conda activate
#移除环境ls
conda deactivate

⑤安装软件

conda install +软件名

⑥小技巧

将conda的常用命令加入全局环境，并改简写，方便使用

vim ~/.bashrc
#添加以下命令至最后一行
alias RNA='conda activate RNA-seq'
#关闭编译器，重启环境即可
source ~/.bashrc

注意：

需要除去默认的defults，并在每个channel末尾加上linux-64/，同时需要末尾添加show_channel_urls: true。

二、PBS文件查看和提交

PBS作业脚本选项(若无-C选项，则每项前面加‘#PBS’)

-a date_time : date_time格式为:[I[[CC]YY]MM]DD]hhmm[.Ss]表示经过date_time 时间后作业才可以运行。

-c interval : 定义作业的检查点间隔，如果机器不支持检查点，则忽略此选项。

-C directive_prefix :在脚本文件中以directive_prefix开头的行解释为qsub 的命令选项。(若无此选项，则默认为'#PBS’)

-e path : 将标准错误信息重定向到 path

-I： 以交互方式运行

-j join: 将标准输出信息与标准错误信息合并到一个文件 join中去。

-k keep : 定义在执行结点上保留标准输出和标准错误信息中的哪个文件。keep为o表示保留前者，e表示后者，oe或eo表示二者都保留，n表示皆不保留。若忽略此选项，二者都不保留。

-l rcsource_list : 定义资源列表。以下为几个常用的资源种类。

cput=N : 请求N秒的CPU时间; N·也可以是 hh:mm:ss 的形式。

mem=N[K|M|G][B|W]: 请求N{kilo/megalgiga } {bytes|words}大小的内存。

nodes=N:ppn=M: 请求N个结点，每个结点M个处理器。

-m mail_options : mail_option为 a: 作业 abort时给用户发信; 为b: 作业开始运行发信；为e: 作业结束运行时发信。若无此选项，默认为a。

-M user_list : 定义有关此作业的 mail发给哪些用户。

-N name: 作业名，限15个字符，首字符为字母，无空格。

-o path: 重定向标准输出到 path。

-p priority: 任务优先级，整数，[-1024，1023]，若无定义则为0.

-q destination : destination有三种形式 : queue , @server,queue@server ，主要是添加运行队列名。

-r y|n: 指明作业是否可运行，y为可运行，n为不可运行。

-S shell: 指明执行运行脚本所用的 shell，须包含全路径。

-u user_list : 定义作业将在运行结点上以哪个用户名来运行。

-v variable_list : 定义export到本作业的环境变量的扩展列表。

-V: 表明qsub命令的所有环境变量都export到此作业。

-W additional_attributes: 作业的其它属性。

-Z : 指明qsub命令提交作业后，不在终端显示作业号。

#较为常用的
#PBS -N default
#PBS -l walltime=336:00:00
#PBS -l nodes=1:ppn=1
#PBS -q cv3
#PBS -j oe
#PBS -o {out}.log
#PBS -V 

三、查看数据质量并进行质控

①fastqc查看数据质量

#文件路径
cd ~/lianxi/rna-seq/fastq/
fastqc *fq.gz
#查看测序量
cat your.fastq.gz | grep -c '^+$'

②获得所有的数据质量后，使用multiqc总的查看

multiqc *.zip

③查看质量值

head SRR851934.1_2.fastq | awk '{if(NR%4==0) printf("%s",$0);}' |  od -A n -t u1 | awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding!";}'

④用trim_galore进行质控

提交循环的脚本run_trim_galore.sh

#循环提交作业

for f in ../clean/*_R1.fq.gz
do
    af=`basename $f`   
    bf=`echo $af | sed s/_R1.fq.gz/_R2.fq.gz/`
    cf=`echo $af | sed s/_R1.fq.gz//g`
    if [ ! -f ${cf}_trimgalore.txt ] #已经清洗过的不需要再提交
    then
        echo $af
        qsub -N default -v pair1=$f,pair2=../$bf,out=$cf do_trimgalore.sh #-v 传参
        sleep 2
    fi
done

运行脚本do_trim_galore.sh

#cutadapt="path"
${trimgalore} -o ./   --length15 --phred33 --stringency 10 --path_to_cutadapt \
#${cutadapt}  测序接头为默认可以不需要
--paired $pair1 $pair2

参数设置：

-o 输出文件夹目录，需要提前建立。

--phred33 选择-phred33或者-phred64，表示测序平台使用的Phred quality score。目前测序基本上都是phred33，根据multiqc结果也能看出来。

--stringency 设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到。

--length：设定输出reads长度阈值，小于设定值会被抛弃。

--paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。

--adapter：输入adapter序列。也可以不输入，Trim Galore会自动寻找可能性最高的平台对应的adapter。自动搜选的平台三个，也可以直接显式输入这三种平台，即--illumina、--nextera和--small_rna。我检查了质控结果，显示illumina通用接头，因此这里不必输入。

--quality：设定Phred quality score阈值，默认为20。

--retain_unpaired：对于双端测序结果，一对reads中，如果一个read达到标准，但是对应的另一个要被抛弃，达到标准的read会被单独保存为一个文件。

--gzip和--dont_gzip：清洗后的数据zip打包或者不打包。

-- trim-n : 移除read一端的reads

④对清洗后的数据再次进行质控

四、基因组数据及注释的下载

这里使用的是UCSC小鼠mm10的基因组数据库以及ensembl的基因组注释文件，但是应该注意一个问题，ensembl数据库的染色体标识为“1”，“2”，“3”等，而UCSC数据库染色体标识为“chr1”，“chr2”， “chr3” 等，二者不匹配，因此需要将ensembl基因注释文件进行修改并清洗。下载方式如下：

wget https://hgdownload.soe.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz
tar -zvxf chromFa.tar.gz
cat *fa>mm10.fa
rm -rf chr*
#获取ensembl基因组数据
wget http://ftp.ensembl.org/pub/release-79/gtf/mus_musculus/Mus_musculus.GRCm38.79.gtf.gz
gunzip -d Mus_musculus.GRCm38.79.gtf.gz

参考：人和小鼠的基因组下载的四种方法：https://www.bilibili.com/read/cv10447213/

数据清洗：

#去除假基因（Pseudogenes），假基因的序列通常与对应的基因相似，但至少是丧失了一部分功能，如基因不能表达或编码的蛋白质没有功能。

VERSION=102
# Changed for v102.
sed 's/^/chr/g' Mus_musculus.GRCm38.102.gtf | sed 's/chrMT/chrM/g' | awk '/^chr[1-9XY]/' | grep 'exon' | grep -E -v 'snoRNA|snRNA|miRNA|rRNA|pseudogene|misc_RNA|_decay|-ps[1-9]|"Rpl[1-9]|"Rps[1-9]' >mm10.ensemblv102.nopsuedo.gtf

五、索引建立以及比对

①使用STAR建立索引以及比对

建立索引

path='~/lianxi/rna-seq'
STAR --runThreadN 4  --runMode genomeGenerate \
--genomeDir $path/index \
--genomeFastaFiles $path/genome/mm10.fa \
--sjdbGTFfile $path/mm10.ensemblv79.nopsuedo.gtf \

针对fastq.gz文件增加--readFilesCommand gunzip -c 参数;针对bzip2文件使用--readFilesCommand bunzip2 -c参数

需要注意“\”后边不能有空格

STAR比对

提交循环的脚本run_star_rsem.sh
mm10index='/public/home/ylsheng_gibh/geodata/mouse/genomeindex/gencodemm10/mm10'
index=${mm10index}
#fastqdir="/public/home/lianxi/rna-seq/fastq/clean"
for f in  *.fq*
do
    root=`basename $f`
    if [[ $root == *.fq]]
    then
        bf=`echo $root | sed -r 's/_R1.fq//g'`
        p1=`echo $root`
        p2=`echo $root | sed 's/_R1.fq/_R2.fq/g'`
        opt='--phred33-quals -p 20 --output-genome-bam --estimate-rspd --star'
    elif [[ $root == *_R1.fq.gz ]]
    then
        bf=`echo $root | sed -r 's/_R1.fq.gz//g'`
        p1=`echo $root`
        p2=`echo $root | sed 's/_R1.fq.gz/_R2.fq.gz/g'`
        opt='--phred33-quals -p 20 --gzipped-read-file --output-genome-bam --estimate-rspd --star'
    fi
    tt=`echo ${bf}.genes.results`
    if ! [ -f $tt ]
    then
        #echo $p1,$p2
        echo PE ... $tt
        qsub -N default -v opt="${opt}",read1=$fastqdir/${p1},read2=$fastqdir/${p2},out=${bf},index=$index do_star_rsem.sh
        sleep 2
        fi
done

运行脚本do_star_rsem.sh
rsem-calculate-expression $opt --paired-end ${read1} ${read2} ${index} ./${out}
rsem-plot-model ${out} ${out}.pdf
rm -rf ${out}.stat

STAR的选项（标注红色为较常用选项）：

--runThreadN：线程数。

--runMode genomeGenerate：构建基因组索引。

--genomeDir：索引目录。（index_dir一定要是存在的文件夹，需提前建好）

--genomeFastaFiles：基因组文件。

--sjdbGTFfile：基因组注释文件。

--twopassMode Basic：使用two-pass模式进行reads比对。简单来说就是先按索引进行第一次比对，而后把第一次比对发现的新剪切位点信息加入到索引中进行第二次比对。

--quantMode TranscriptomeSAM GeneCounts：将reads比对至转录本序列。

--alignIntronMin：最短的内含子长度。（根据GTF文件计算）

--alignIntronMax：最长的内含子长度。（根据GTF文件计算）

--outSAMtype BAM SortedByCoordinate：输出BAM文件并进行排序。

--sjdbOverhang：reads长度减1。

--outSAMattrRGline：ID代表样本ID，SM代表样本名称，PL为测序平台。在使用GATK进行SNP Calling时同一SM的样本可以合并在一起。

--outFilterMismatchNmax：比对时允许的最大错配数。

--outSJfilterReads Unique：对于跨越剪切位点的reads（junction reads)，只考虑跨越唯一剪切位点的reads。

--outSAMmultNmax：每条reads输出比对结果的数量。

--outFileNamePrefix：输出文件前缀。

--outSAMmapqUnique 60：将uniquely mapping reads的MAPQ值调整为60，满足下游使用GATK进行分析的需要。

--readFilesCommand：对FASTQ文件进行操作。

--readFilesIn：输入FASTQ文件的路径。

获得表达矩阵

使用rsem输入转录组的bam文件进行差异表达定量

②使用rsem-STAR进行差异定量

建立索引

使用rsem调用STAR进行索引的建立

path=~/lianxi/rna-seq
rsem-prepare-reference -gtf $path/genome/mm10.ensemblv79.nopsuedo.gtf \
--STAR $path/genome/mm10.fa \
$path/index/mm10/mm10

得出报告文件是STAR默认的索引建立，有以下小问题：

sjdbOverhang选项是100，而illumina测序序列长度是150，经过搜索得出结论是在大多数情况下，通用值100和理想值一样有效。

差异分析

报错：这里两个mm10，第一个是文件夹，第二个是基因名称，只提交了一个被报错。

或者使用shell脚本调用rsem-STAR流程分析

将比对结果分析图绘制pdf，需要安装Rscript

cd ~/lianxi/rna-seq/mapping
ls *.log|while read id;do (rsem-plot-model ${id%%.*} ${id%%.*}.pdf );done
rm -rf *.stat

可以在上一步直接加在末尾就不用重写循环。

最后结果：

index结果：

rsem-STAR结果

  

mm10index='/public/home/ylsheng_gibh/geodata/mouse/genomeindex/gencodemm10/mm10'
index=${mm10index}
#fastqdir="/public/home/lianxi/rna-seq/fastq/clean"
for f in  *.fq*
do
        root=`basename $f`
        if [[ $root == *.fq]]
        then
        bf=`echo $root | sed -r 's/.fq//g'`
        opt='--phred33-quals -p 20 --output-genome-bam --estimate-rspd --star'
        fi
tt=`echo ${bf%%.*}_trimmed.genes.results`
        if ! [ -f $tt ]
        then
                echo $p1,$p2
                echo PE ... $tt
                qsub -N default -v opt="${opt}",read=/public/home/ylsheng_gibh/geodata/mouse/clean/$f,out=${bf},index=$index do_star_rsem.sh
                sleep 2
        fi
done