RNA-seq全流程基础版之上游操作

全流程介绍

1. 公共数据库下载数据
2. 转换数据库格式成为fasta.gz
3. fastqc或multiqc质控，trim-galore去接头和低质量序列
4. hisat2比对到基因组获取sam文件
5. samtools将sam文件转化为bam文件
6. subread进行count计数得到表达矩阵

软件安装(linux系统)
安装miniconda:wget -c https://mirrors.tuna.tsinghua.edu.cn/anaconda/miniconda/Miniconda3-latest-Linux-x86_64.sh 清华源经常崩溃可自行更换镜像源。
安装后
source .bashrc
配置镜像

conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/main/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/pkgs/free/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/conda-forge/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/msys2/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/bioconda/
conda config --add channels https://mirrors.ustc.edu.cn/anaconda/cloud/menpo/
conda config --set show_channel_urls yes

miniconda安装后下列软件均可用conda安装也可以上官网找到对应版本链接wget安装
安装sra-tools、fastqc、multiqc、trim-galore、hisat2、samtools、subread。

conda install -y sra-tools # 出现三个done,即表示安装成功，否则重新提交命令，下同。
conda install -y fastqc
conda install -y fastqc
conda install -y trim-galore
conda install -y samtools
conda install -y subread

安装比对软件；hisat2并建立索引（常见物种有官方索引）

conda install -y hisat2

wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg38.tar.gz

**下载数据 **
sra-tools下载单个数据

prefetch SRR号

sra-tools批量下载数据：

nohup prefetch --option-file SRR_Acc_List.txt &

将sra文件转换成fasta.gz文件：fasta-dump或fasterq
常见参数：
-O|–outdir ：指定输出目录 ./当前文件夹（默认）
-e|–threads：指定线程数，默认为6
-p|–progress ：显示进度
-s|–split-spot ：双端测序结果存储在1个文件
-S|–split-files：双端测序结果存储在2个文件
-3|–split-3： 双端测序结果存储在3个文件 （第三个文件存放没有配对的read）
批量用法，举例三种

for id in *sra; do pfastq-dump --threads 4 ./$id --gzip; done

for id in *sra; do pfastq-dump --threads 4 ./$id --split-3 --gzip; done

for id in *sra
do
echo $id
fastq-dump --split-3 --gzip $id
done

质控
fastqc（生成两个文件* .html和*fastqc.gz），multiqc

ls SRR*.gz | while read id;do (nohup fastqc $id &);done

for id in *.gz
do 
echo $i
nohup fastqc $id -t 4 -o ./ &
done

multiqc：使用依赖python3（创建python=3的rnaseq2环境）

multiqc ./#对当前目录所有文件进行multiqc

去接头和低质量序列
trim-galore常用参数。
–quality：设定Phred quality score阈值，默认为20。分析时可改成25，稍微严格一些。
–phred33：选择-phred33或者-phred64，表示测序平台使用的Phred quality score。一般phred33对应(Sanger/Illumina 1.9+ encoding)，-phred64对应(Illumina 1.5 encoding)
–stringency：设定可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。
–length：设定输出reads长度阈值，小于设定值会被抛弃。
–paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准。
–gzip和–dont_gzip：清洗后的数据zip打包或者不打包。
–output_dir：输入目录。需要提前建立目录。

for（ i=71;i<=80；i++）in SRR144284$i.fastq.gz
do 
echo $i
nohup trim_galore -q 25 --length 36 --phred33 --paired --gzip -o ./clean/ &
done

ls *_1.fastq.gz >1
ls *_2.fastq.gz >2
paste 1 2 > config
cat config|while read id
do
        arr=($id)
        fq1=${arr[0]}
        fq2=${arr[1]}
nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 3 --gzip --paired -o ./clean/ $fq1 $fq2 &
done

for ((i=71;i<=80;i++));trim_galore -q 25 --length 36 --phred33 --paired --gzip -o ./clean/ &;done

比对到参考基因组（获取sam文件）
hisat2（有官方的索引文件）
建立基因组索引：

hisat2-build –p 4 genome.fa genome

建立基因组，转录组，SNP索引：

hisat2-build -p 8 genome.fa --snp genome.snp --ss genome.ss --exon genome.exon genome_snp_tran_index

hisat2用法：hisat2 -x 索引文件的前缀 -1 第一个文件 -2 第二个文件 -S 指定输出路径
hisat2常用参数。
-x：参考基因组索引文件的前缀。
-1：双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2：双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
-U：单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
–sra-acc：输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
-S：指定输出的SAM文件。

 for ((i=71;i<=80;i++));
do nohup hisat2 -t -x ./reference/grch38/genome \
-1 ./rna_seq/data/SRR144284${i}.sra_1.fastq.gz \
-2 ./rna_seq/data/SRR144284${i}.sra_2.fastq.gz \
-S SRR144284${i}.sam &;done

将sam转化成bam并排序（bam是sam文件的二进制形式，更加节约空间）

for ((i=56;i<=62;i++));do samtools sort SRR144284${i}.bam -o SRR144284${i}_sorted.bam;done

hisat2比对的同时转换成bam文件

index=./reference/grch38/genome
ls *_1_val_1.fq.gz|while read id
do
id=${id/_1_val_1.fq.gz/}
hisat2 -p 1 -x ${index} \
-1 ${id}_1_val_1.fq.gz \
-2 ${id}_1_val_1.fq.gz \
2>${id}.hisat2.log| \
samtools sort -@ 2 -o ${id}.hisat2.sorted.bam
done

count计数
subread（含featureCounts）
subread常用参数
-input file：输入的bam/sam文件，支持多个文件输入
-a < string > ：参考gtf文件名，支持Gzipped文件格式
-F ：参考文件的格式，一般为GTF/SAF，C语言版本默认的格式为GTF格式
-A：提供一个逗号分割为两列的文件，一列为gtf中的染色体名，另一列为read中对应的染色体名，用于将gtf和read中的名称进行统一匹配，注意该文件提交时不需要列名
-g：gene_id表示meta-feature名称为gene_id(ensembl名称)
-M：如果设置-M，多重map的read将会被统计到
-o < string >： 输出文件的名字
-O：允许多重比对，即当一个read比对到多个feature或多个metafeature的时候，这条read会被统计多次
-p：双端测序时使用
-T ：线程数目，1~32

运行featureCounts计数

featureCounts -T 6 -p -t exon -g gene_id -a ./gtf/gencode.v25.annotation.gtf.gz -o ./all.id.txt *sort.bam

简化结果，去除特征列，只保留基因计数的列

cat all.id.txt | cut -f1,7- > counts.txt