学习资料：B站 jimmy大神的视频
一、质控命令
用到的fastqc是北鲲云上自带的，multiqc是在conda的py3.7环境里自己装的。
针对昨天转成fastq.gz的四个数据进行质控：
建立rnaseq-qc.slurm文件。对比以前在学校服务器上编的sh文件，其实没有什么太多不同。只不过就是要调用自带软件需要输入module add XXX 罢了，别的都一样。
#!/bin/bash
#SBATCH --output=rnaseq-qc.out
#SBATCH --error=rnaseq-qc.err
#SBATCH --mail-type=end
#SBATCH --mail-user=zmeraner@126.com
module add Anaconda3/2020.02
source activate
conda activate py3.7 #激活3.7环境
module add FastQC/0.11.9-Java-11 #调用fastqc软件
project=~/rnaseq #指定项目位置
mkdir -p $project/qc/raw
ls $project/rawdata/*.gz | xargs fastqc -t 10 -O $project/qc/raw #将rawdata中的数据逐一进行质控
multiqc $project/qc/raw -o $project/qc/raw #将raw文件夹里的质控文件整合
fastqc生成下面8个文件，multiqc生成上面2个文件。
image.png
二、质控文件解读
fastqc文件主要看总测序量reads数(一般表达谱20M左右就够，后面详述)，还有质量（Per base sequence quality最好>30）
GC含量等，如果为红色或橙黄色表明有问题，绿色表示这个指标没问题。比如这里面有一些项目似乎有点儿问题，应该查看一下。
image.png
per base sequence content这个参数理论上来说，A和T应该相等，C和G应该相等，且相互平行，但可能测序一开始不稳定，刚开始测序仪状态不稳定，所以会出现这种波动的情况，需要进行修剪。
看一下multiqc里面，四个样本都是红色的，都需要对前面1-15个bp进行修剪。
image.png
multiqc的文件可以在一个文件里看所有样本的质量数据。
image.png
duplicate这个参数4个样本都是橙黄色，说明需要去重。
image.png
三、测序量
测序量
前面在质控数据里有个Total Sequences，这个是总共有多少reads。
对于单端测序的PE50，total sequnceses是20M的话，总测序量应该是20M50,也就是1G。
现在送测序多为PE150，一般测序reads也要达到20M，那么测序量就是20M300，总共就是6G。所以一般自己建库好了，送测序就测6G，一个G就几十块钱，真不贵啊。回来自己看看total sequences，到了20M左右，就表明测序量满意。想知道到底多少测序量是够的呢，找了些资料看看：
如果是就看看表达谱那么每个样本可能只需要5–25 M reads
如果要得到很全面的表达信息或可变剪接，每个样本需要30–60 M reads
如果要深度研究或者新物种可能要每个样本100–200 M reads
小RNA或者靶向RNA测序每个样本需要1–5 M reads
测序深度
人类基因组有3000Mnt，其中大约1/30被用于蛋白质编码基因。要测序一遍RNA大约需要100M 数据量。那么测序量÷100M，大约就是测序深度。
比如测了6G，6G÷100M=60X，意思是RNA被测了60遍。但是一般DNA常常要说测序深度，而RNA，因为不同基因的表达量不同，而且长度不同，所以不是所有的基因都被测了这么多啊。
四、trim
在py3.7环境下安装trim-galore conda install trim-galore
建立文件rnaseq_trim.slurm
输入下面命令：用trim-galore进行修剪，再进行一遍质控。
#!/bin/bash
#SBATCH --output=rnaseq_trim.out
#SBATCH --error=rnaseq_trim.err
#SBATCH --mail-type=end
#SBATCH --mail-user=zmeraner@126.com
module add Anaconda3/2020.02
source activate
conda activate py3.7 #激活3.7环境
project=~/rnaseq #项目文件夹
#step3：trim & qc
mkdir $project/clean
output_dir="$project/clean"
ls $project/rawdata | grep "fastq.gz" > config_file
cat config_file | while read id
do
trim_galore --clip_R1 15 --phred33 --quality 20 --length 20 --stringency 3 -o $output_dir $project/rawdata/$id
done
module add FastQC/0.11.9-Java-11
mkdir -p $project/qc/clean
ls $project/clean/*.gz | xargs fastqc -t 10 -O $project/qc/clean
multiqc $project/qc/clean -o $project/qc/clean
--clip_R1 15这个参数是把5`端15个bp去掉，现在长度应该是35bp了。
image.png
生成的文件带有trimmed标记，也是fq.gz格式
image.png
网友评论