生信技能树学习笔记
Raw data
背景
先了解
- 测序长度
- 单端/双端？
- 测序对象 mRNA？lncRNA?
fastq数据格式
Raw data 或 Raw reads 结果以FASTQ文件格式存储
结果每四行一显示
- 第一行 @开头，随后为illumina测序识别符合描述文字
- 第二行 碱基序列
- 第三行 开头
- 第四行 对应序列的测序质量的ASCII码 Base calling，Q值越大精度越高，ASCII数值减33得到Q值
质控——fastqc
常用参数
- -o --outdir 设置输出目录
- -t --threads 同时处理几个样本
nohup fastqc -t 6 -o ./ SRR*.fastq.gz >qc.log &将html文件下载到本地，解读 QC 报告
- Basic Statistics——GC 30-60%
- Per base sequence quality——横坐标为碱基位置，纵坐标为Q值
- Per tile sequence quality——好的结果应该是全蓝色
- Per sequence quality scores——横坐标为平均Q值，纵坐标为每个Q值对应的reads数（应该呈现拖尾的分布）
- Per base sequence content——每个碱基位置上ATCG含量的分布图，AT和GC应分别相等，呈水平线，开头允许少许抖动
- Per sequence GC content——横坐标为平均GC含量，纵坐标为每个GC含量对应的序列数量，蓝色为理论值，红色为测量值，二者越接近越好
- Per base N content——N含量分布
- Sequence Length Distribution——长度分布
- Sequence Duplication levels——序列的重复度
- Overrepresented sequences——转录组中某个
- Adapter Content——接头含量，表示序列中两端adapter的情况
使用MultiQC整合FastQC结果
代码语言：{ r setup, include = FALSE}复制multiqc *.zip -o ./ 数据过滤
如何判断数据需要过滤？
质量控制标准
- 去除含接头的reads
- 过滤去除低质量值数据，确保数据质量
- 去除含有N（无法确定碱基信息）的比例大于5%（根据实际情况）的reads
数据过滤方式一：trim_galore
常用参数
- -q --quality 切除质量得分低于设置值的序列，默认值20
- --length 长度小于设定值的reads将被丢弃
- --max_n 去除含有碱基数大于N的序列
- --stringency 限定最少与adaptor序列重叠的碱基数
## 单个样本
trim_galore -q 20 --length 20 --max_n 3 --fastqc --paired -o ./ ../../rawdata/SRR1039510_1.fastq.gz ../../rawdata/SRR1039510_2.fastq.gz## 多个样本
## 生成并保存ID，$NF表示最后一列
ls $HOME/project/Human-16-Asthma-Trans/data/rawdata/*_1.fastq.gz | awk -F'/' '{print $NF}' | cut -d'_' -f1 >ID
vim trim.sh
## 循环 .sh内容
rawdata=$HOME/project/Human-16-Asthma-Trans/data/rawdata
cleandata=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
cat ID | while read id
do
trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata} ${rawdata}/${id}_1.fastq.gz ${rawdata}/${id}_2.fastq.gz
done
## 后台运行
nohup sh trim.sh > trim.sh.log &
## 整合新的multiqc
multiqc *.zip数据过滤方式二：fastp
速度比 trim_galore 快
常用参数
- -i, -I 后接需要过滤的fastq文件
- -o,-O 后接过滤玩输出的fastq文件名 【注意大小写和reads1/2前后对应】
- -n --n_base_limit 限制N个数
- -q 设置碱基质量阈值，默认阈值为15
- -l 小写的L 设置 read 的最小长度，默认是15，长度＜15 的 read 将被丢弃
# 定义文件夹：vim fastp.sh
cleandata=$HOME/project/Human-16-Asthma-Trans/data/cleandata/fastp/
rawdata=$HOME/project/Human-16-Asthma-Trans/data/rawdata/
cat ../trim_galore/ID | while read id
do
fastp -l 20 -q 20 --compression=6
-i ${rawdata}/${id}_1.fastq.gz
-I ${rawdata}/${id}_2.fastq.gz
-o ${cleandata}/${id}_clean_1.fq.gz
-O ${cleandata}/${id}_clean_2.fq.gz
-R ${cleandata}/${id}
-h ${cleandata}/${id}.fastp.html
-j ${cleandata}/${id}.fastp.json
done
## 注意反斜线后面一定不能有空格
# 运行fastp脚本
nohup sh fastp.sh >fastp.log &过滤数据比较
代码语言：{ r setup, include = FALSE}复制# 进入过滤目录
cd $HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
# 原始数据
zcat $rawdata/SRR1039510_1.fastq.gz | paste - - - - > raw.txt
# 过滤后的数据
zcat SRR1039510_1_val_1.fq.gz |paste - - - - > trim.txt
awk '(length($4)<63){print$1}' trim.txt > Seq.ID
head -n 100 Seq.ID > ID100
grep -w -f ID100 trim.txt | awk '{print$1,$4}' > trim.sm
grep -w -f ID100 raw.txt | awk '{print$1,$4}' > raw.sm
paste raw.sm trim.sm | awk '{print$2,$4}' | tr ' ' 'n' |less -S附录
前台转后台
代码语言：{ r setup, include = FALSE }复制## 运行时按Ctrl Z
jobs ## 可以看到 stopped，并且可以看到 jobs 的序列号，但只能在当前窗口
ps fx ## 可以看到PID编号，使用ps fx在另一个窗口也可以看见
bg %1 ## 百分号后面的是jobs的序列号
jobs ## 此时进入running状态，并且可在后台显示
kill -9 %2 ## 强制性杀除第2个进程
fg %1 ## 此时后台转到前台检查脚本内容
使用echo打印命令，进行检查
代码语言：{ r setup, include = FALSE}复制cat ID | while read id
do
echo "trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata} ${rawdata}/${id}_1.fastq.gz ${rawdata}/${id}_2.fastq.gz"
done掐头去尾获得样本ID
代码语言：{ r setup, include = FALSE}复制ls $rawdata/*_1.fastq.gz | while read id
do
name=${id##*/}
name=${name%_*}
trim_galore -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ${cleandata} ${rawdata}/${name}_1.fastq.gz ${rawdata}/${name}_2.fastq.gz
done
