RNA-seq的上游分析主要基于Linux平台，根据给定的fastq文件，经过系列生信软件处理最终得到表达矩阵的过程。原始数据文件一般可以从文献中获取（可能是sra格式，需要转换成fastq，之后会专门整理下如何获取）。我之前练习的是老师直接给我的fastq文件(猴子的数据)，还要补充的是双末端测序同一样本会有两个文件的。就从这开始吧~
0、环境配置所需软件
conda activate rna-seq
conda install fastqc -y
conda install multiqc -y
conda install trimmomatic -y
conda install hisat2 -y
conda install trim-galore -y
conda install samtools -y
conda install -c bioconda subread
- 关于conda的说明，见之前的整理
1、质控
目的：检测fastq测序质量如何，不好的话要进行一定的处理才能用以分析
涉及软件：fastqc、multiqc
ls *gz | xargs fastqc -o /home/monkey4/myfastqc -t 10
- 注意此时所在位置要在fastq文件路径下
-
xargs命令实现了将ls的结果作为fastqc的参数输入 - fastqc的
-o选项交代结果储存路径；-t选项交代线程数 - 关于qc的结果，每个fastq文件都会有一个对应的“检查报告”，也可以将所有报告进行合并。
cd /home/monkey4/myfastqc
multiqc ./
# multiqc会进行自动识别合并
- 关于报告结果解读，网上有很多总结的介绍。之后自己可能也会整理下。目前的测序结果都蛮不错的，比如我做的那个。
- 如果的确测序结果质量不行，可以进行下面的处理。（我自己没有实操过..仅记录下）
#单独处理
trim_glaore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 --paired -o 输出路径 路径/xxx_1.fastq.gz 路径/xxx_2.fastq.gz
#批量处理
ls 路径/*_1.fastq.gz >1
ls 路径/*_2.fastq.gz >2
paste 1 2 #查看
paste 1 2 > config
cat >qc.sh
conda activate rna-seq
bin_trim_glore=trim_glaore
dir='/home/li/monkey4/test/clean'
cat $1 |while read id #注意这里的$1 表示第一个参数，即bash qc.sh config 中的config
do
arr=($id)
fq1=${arr[0]}
fq2=${arr[1]}
nohup $bin_trim_glore -q 25 --phred33 --length 36 -e 0.1 --stringency 3 -paired -o $dir $fq1 $fq2 &
done
source deactivate
ctrl+c #终止输入编辑
bash qc.sh config
2、比对
- 目的：将fastq文件的所有read短序列比对到参考基因组上，即read本来属于哪条基因上的。
- 涉及软件：hisat2（还有同类型软件bowtie2等）、samtools
2.1、构建索引
构建索引的原因是因为参考基因组文件太太太大了，把一个几百bp的read比对到基因组上，难度可想而知。而索引文件就类似图书的目录，按图索骥就方便很多。
hisat2-build -p 10 基因组文件路径/xxxx.fa 索引储存路径/xxxx
#索引文件命名不需要加后缀，取个名即可。软件会自己加的
hisat2-build -p 4 /home/li/genome/monkey.fa /home/li/monkey4/hisat2-index/Macaca
#这个是我的实操代码
2.2、进行比对
#单个比对
hisat2 -t -p 4 -x 索引文件路径/索引文件名 -1 XXX1.fq.gz -2 XXX2.fq.gz -S 储存路径/XXX.sam
# 针对双末端测序，可以同时比对，进行结果合并
hisat2 -t -p 4 -x /home/li/monkey4/rnaseq_test/index/hisat2_index/Macaca -1 M5Acortex_1.fq.gz -2 M5Acortex_2.fq.gz -S /home/li/monkey4/rnaseq_test/M5A.sam
#实操代码
#批量比对
ls *gz |cut -d"_" -f 1 | sort -u |while read id ; do hisat2 -p 10 -x 索引文件路径/索引文件名 -1 ${id}_1.fq.gz -2 ${id}_2.fq.gz -S 输出路径/${id}.hisat.sam ; done
ls *gz |cut -d"_" -f 1 | sort -u |while read id ; do hisat2 -p 10 -x /home/li/monkey4/rnaseq_test/index/hisat2_index/Macaca -1 ${id}_1.fq.gz -2 ${id}_2.fq.gz -S /home/li/monkey4/rnaseq_test/sam/${id}.sam 1>$id.download.log 2>&1 ; done
分析得到的sam文件记录一个样本的所有read比对信息。
2.3、sam转bam
samtools view -S XXX.sam -@ 10 -b > XXX.bam
ls *.sam | cut -d"." -f 1 | sort -u | while read id ;do samtools view -S ${id}.hisat.sam -@ 10 -b > ${id}.bam ; done
而bam格式中的b是binary的意思，是sam格式的二进制表示方式，主要就是为了节省内存，提高后面分析速度。
3、计数
- 目的：统计每个基因比对上了多少个read
- 涉及软件：subread
#单个计数
featureCounts -T 10 -p -t exon -g gene_id -a 基因组注释文件路径/XXX.gtf -o 结果文件名 bam文件路径/XXX.bam
featureCounts -T 10 -p -t exon -g gene_id -a /home/li/genome/Macaca_fascicularis.Macaca_fascicularis_5.0.97.gtf -o counts.M2Ccortex.txt /home/li/monkey4/test/aligned/M2Ccortex.bam
#批量处理
ls *tex.bam |cut -d"." -f 1 |sort -u |while read id ; do featureCounts -T 10 -p -t exon -g gene_id -a /home/li/genome/Macaca_fascicularis.Macaca_fascicularis_5.0.97.gtf -o counts.${id}.txt /home/li/monkey4/test/aligned/${id}.bam; done
#将批量计数结果转成CSV格式
ls *.txt |cut -d"." -f 1,2 |sort -u | while read id ; do cut -f 1,7-12 ${id}.txt | grep -v "#" > countmatrix.${id}.csv; done
#合并所有样本的CSV文件
paste *.csv | awk '{printf $1 "\t" ; for (i=2; i<=NF; i=i+2) printf $i "\t" ; print $i}' > merge.count.csv
以上就是一个极简版的RNA-seq的上游分析过程，最终获得了基因表达矩阵，用以后期的数据分析。文中提到一些数据文件及格式还是值得后期慢慢梳理的。暂时先这样吧。上述代码主要参考自生信技能树的学习，站在巨人肩膀学生信，嘻嘻~
