exit1是退出的意思
忽视输入的报错，重定向在日志文件，不是报错【error】
image.png
image.png
转录组分析总流程
image.png
数据质控
背景知识
-
数据量的统计方式image.png
sra转换成fastq
image.png
# 定义文件夹
fqdir= ~/project/Human_16-Asthma-Trans/data/rawdata/sra/fastq
# 单个转换
fastq-dump --gzip --split-2 -X 25000 -O ${fqdir} SRR1039510.sra
#fasterq-dump --split-files SRR11180057.sra
##-O是输出到命令
标准输出格式；为什么要拆分呢？因为后续的软件不支持sra格式
image.png
# 批量转换
ls /trainee2/Nov10/project/Human_16-Asthma-Trans/data/rawdata/sra/*sra|while read id
do
echo "fastq-dump --gzip --split-e -X 25000 -O ${fqdir} ${id}"
done >sra2fq.sh
# 提交后台运行命令
nohup sh sra2fq.sh >sra2fq.log &
质控
image.png
image.png
image.png
- But what does this quality score mean?
The quality score for each sequence is a string of characters, one for each base of the nucleic sequence, used to characterize the probability of mis-identification of each base.- The score is encoded using the ASCII character table :
image.png
image.png
image.png
image.png
image.png
image.png
显示第一行
质控软件
image.png
# 激活conda环境
conda activate rna
# 连接数据到自己的文件夹
ln -s /teach/data/airway/fastq_raw25000/*gz .
# 使用FastQC软件对单个fastq文件进行质量评估，结果输出到qc/文件夹下
qcdir=~/project/Human_16-Asthma-Trans/data/rawdata/qc
fqdir=~/project/Human_16-Asthma-Trans/data/rawdata/fastq
fastqc -t 3 -o $qcdir $fqdir/SRR1039510_1.fastq.gz
image.png
生成的文件通过Xftp拖出
双击打开html后文件如下
如何看质控报告
绿色的勾勾表示达到了标准，xx表示没有达到标准
-
Basic Statistics
image.png
image.png -
数据量的统计方式
通常说数据量多少是指的测了多少个碱基，而不是说这个文件的大小
3.9G是测序的数据量，2.2G是文件的大小
-
Per base sequence quality
每一个位置的的碱基都有一个对应的箱图
-
Per tile sequence quality
image.png
-
Per Sequence Quality Scores
image.png
-
Per Sequen GC
image.png
# 多个数据质控
fastqc -t 2 -o $qcdir $fqdir/SRR*.fastq.gz
# 使用MultiQc整合FastQC结果
multiqc *.zip
image.png
数据过滤
-
过滤条件
image.png
1、adapter是一段短的序列已知的核酸链，用于链接序列未知的目标测序片段。
2、barcode，也称为index，是一段很短的寡居核酸链，用于在多个样品混合测序时，标记不同的样品。
3、insert是用于测序的目标片段，因为是包括在两个adapter之间，所以被称为“插入”片段。
一个常见测序片段类似与adapter--barcode--insert--adapter。测序开始时前几个碱基无法测得，第一个adapter在数据输出时被去除；由于测序仪读长限制，第二个adapter通常无法测得。所以，经常得到类似 barcode--部分insert的read。最后，把barcode去除，只保留测度insert的片段，这个操作的术语是demultiplexing。但是有时候测序时会测穿，也就是说会得到barcode--insert的read--部分adapter，那么这里就包含了接头了，这里的接头也就是大家经常说去接头要去除的部分。
image.png
image.png
-
trim_galore过滤
image.png
# 定义文件夹
rawdata=/teach/project/Human-16-Asthma-Trans/data/rawdata/fastq
cleandata=/teach/project/Human-16-Asthma-Trans/data/cleandata/trim_galore
# 单个样本
trim_galore --phred33 -q 30 --length 30 --stringency 3 --fastqc --paired --max_n 3 -o $cleandata $rawdata/SRR1039510_1.fastq.gz $rawdata/SRR1039510_2.fastq.gz
##length是小于30的不要 max——n指的是含有未知的大于3就要去掉了
image.png
image.png
# 多个
cat /teach/project/Human-16-Asthma-Trans/data/rawdata/sra/sampleId.txt | while read id
do
echo "trim_galore --phred33 -q 20 --length 36 --stringency 3 --fastqc --paired --max_n 3 -o ${cleandata} ${rawdata}/${id}_1.fastq.gz ${rawdata}/${id}_2.fastq.gz"
done >trim_galore.sh
nohup sh trim_galore.sh >trim_galore.log &
# 使用MultiQc整合FastQC结果
multiqc *.zip
-
fastp过滤
image.png
# 定义文件夹
cleandata=/teach/project/Human-16-Asthma-Trans/data/cleandata/fastp
# 单个样本
fastp -i $rawdata/SRR1039510_1.fastq.gz -I $rawdata/SRR1039510_2.fastq.gz \
-o $cleandata/SRR1039510_1.fastp.fq.gz -O $cleandata/SRR1039510_2.fastp.fq.gz \
-l 36 -q 20 --compression=6 -R $cleandata/SRR1039510 \
-h $cleandata/SRR1039510.fastp.html -j $cleandata/SRR1039510.fastp.json
# 多个样本
cat /teach/project/Human-16-Asthma-Trans/data/rawdata/sra/sampleId.txt | while read id
do
echo "fastp -i ${rawdata}/${id}_1.fastq.gz -I ${rawdata}/${id}_2.fastq.gz -o ${cleandata}/${id}_1.fastp.fq.gz -O ${cleandata}/${id}_2.fastp.fq.gz -l 36 -q 20 --compression=6 -R ${cleandata}/${id} -h ${cleandata}/${id}.fastp.html -j ${cleandata}/${id}.fastp.json 1>$cleandata/${id}.fastp.log 2>&1"
done >fastp.sh
# 运行fastp脚本
nohup sh fastp.sh >fastp.log &
-
数据过滤前后的比较
# 进入过滤目录
cd /teach/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 > ID
head -n 100 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
