优秀学员笔记：转录组课堂笔记Day1~2

Day0 背景介绍

Markdown使用教程：

• /

Day1 准备工作

一、工作目录创建

1.1 目录分类

通常存放在home目录，创建子目录以分门别类数据，常见情况如下：

• mkdir project 存放数据分析及其结果
• mkdir project_backup 存放备份数据
• mkdir tools 存放小代码工具
• mkdir pipeline 存放分析流程
• mkdir database 存放背景数据
• mkdir biosoft 存放软件包

注：善用 “_” 进行文件命名

1.2 创建流程

# 进入到个人home目录
cd ~
pwd

# 1.建立数据库目录
mkdir database 

# 2.建立并进入项目目录
mkdir project
cd project

# 3.建立并进入子项目目录
# 注意项目命名习惯：物种-样本数-疾病-分析流程/scRNA/spatial/srna/mirna/circrna
mkdir Human-16-Asthma-Trans 
cd Human-16-Asthma-Trans

# 4.建立数据存放目录(同时建立data及其下rawdata和cleandata目录)
mkdir -p data/rawdata data/cleandata 
# 建立比对目录
mkdir -p Mapping
# 建立定量目录
mkdir -p Expression
# 创建差异分析目录
mkdir Diff_Analysis

# 5.使用tree命令查看目录结构
# 默认展开所有层，使用-L参数控制展开的目录层级（1层）
tree ./ -L 1 

二、准备分析数据

2.1 链接数据

# 激活conda环境
conda activate rna

# 连接数据到自己的文件夹
cd $HOME/project/Human-16-Asthma-Trans/data/rawdata
pwd

# 链接（-s软链接，类似创建快捷方式）
ln -s /home/t_rna/data/airway/fastq_raw25000/*gz  ./

2.2 FASTQ文件介绍

双端测序：一般一个样本对应两个fq文件，gz是压缩后缀，如SRR1039510对应

• read1：SRR1039510_1.fastq.gz
• read2：SRR1039510_2.fastq.gz

FASTQ格式文件中每个Read由四行描述：

• 第一行以“@”开头，随后为llumina测序识别符(Sequence ldentifiers)和描述文字(可选)
• 第二行是碱基序列
• 第三行以“+”开头，随后为Ilumina测序识别符(可选)
• 第四行是对应序列的测序质量的ASCII码(碱基质量值)

碱基质量值(Quality Score或Q-score)是碱基识别(Base Calling)出错概率的整数映射，碱基质量值越高表明碱基识别越可靠，准确度越高。通常使用的碱基质量值Q公式为：

P为碱基识别出错的概率，对应关系如下：

2.3 随堂练习

1. 统计SRR1039510_1.fastq.gz文件中共有多少条reads

# 统计包含@SRR的行数
zless -SN SRR1039510_1.fastq.gz | grep "@SRR" -c
zless -SN SRR1039510_1.fastq.gz | grep '^@SRR' |wc -l

# 将文件内容按行并排合并，每行4列，使用制表符（\t）分隔，统计行数
zless -SN SRR1039510_1.fastq.gz | paste - - - - |wc -l

# 统计总行数，FASTQ 格式每4行表示一个序列条目，总行数除以4得到序列条目数。
zless -SN SRR1039510_1.fastq.gz |wc -l | awk '{print $0/4}'
zless -SN SRR1039510_1.fastq.gz |awk '{ if(NR%4==2) {print} }' |wc -l

# sed '1~4p' 表示从第1行开始，每隔4行打印一次
zless SRR1039510_1.fastq.gz | sed -n '1~4p' |wc -l

1. 输出SRR1039510 1.fastq.gz文件中所有的序列ID(第一行)

zless SRR1039510_1.fastq.gz | grep '^@SRR' |less -S
zless SRR1039510_1.fastq.gz | paste - - - - |cut -f 1 |less -S
zless SRR1039510_1.fastq.gz | awk '{if(NR%4==1){print}}' |less -S

1. 输出SRR1039510 1.fastq.gz文件中所有的序列(第二行)，并统计碱基数

# 取第二行，删除末尾换行符后统计
zless -S SRR1039510_1.fastq.gz | paste - - - - | cut -f 2 | tr -d '\n' | wc -m

# 取出所有类型的碱基“A/T/C/G/N”排列成行，统计行数
zless -S SRR1039510_1.fastq.gz | paste - - - - | cut -f 2 | grep -o [ATCGN] | wc -l

1. 统计SRR1039510 1.fastq.gz文件中的重复序列，从高到低排序

zless -S SRR1039510_1.fastq.gz | sed -n '2~4p' | sort | uniq -c | sort -n k1 | head

三、数据质量评估

3.1 FastQC

[FastQC](Babraham Bioinformatics - FastQC A Quality Control tool for High Throughput Sequence Data[1])可以对Fastq格式的原始数据进行质量统计，评估测序结果，为下一步修剪过滤提供参考。

# 使用FastQC软件对单个fastq文件进行质量评估，结果输出到qc/文件夹下
# 脚本后台运行(nohup &，nohup表示no hang up退出终端不影响，&表示后台运行)
nohup fastqc -t 6 -o ./ SRR*.fastq.gz 1>qc.log 2>&1 &
# 1 表示标准输出，2表示错误输出
# 2>&1 表示将写入2中的日志重定向到1中
# 1>qc.log 2>&1 表示将1和2的任务日志都输出到qc.log中

3.2 如何查看后台任务进度

htop 与 ps fx 可査看任务进程PID，已结束的任务不会显示

①htop调出服务器资源管理器

②ps fxww查看任务投递窗口

③jobs仅查看当前任务

• Running表示任务正在运行
• Done表示任务运行成功并结束
• Exit表示任务失败退出并结束

④查看任务日志

Day2 数据质控与过滤

一、报告解读

1.1 基础统计信息

• Filename：原始分析文件名
• File type：文件类型 base calls 表示碱基识别，带有 fq 序列信息文件
• Encoding：ASCII编码的质量值系统
• Total Sequences：测序序列的总read数
• Total Bases：以碱基(对)数或者测序read数统计测序数量的大小（与计算机存储单位区分）

• Sequence Length：提供最长和最短read的length，如多个read同样长则只显示一个数
• %GC：GC含量，通常在40~60

1.2 序列质量值 fq文件中每一个位置上的所有read的碱基质量值的箱线图。横坐标表示碱基的位置，纵坐标表示碱基的QC值。 在测序过程中，随着测序读长（read length）的增加，碱基质量值下降的现象主要有以下几个原因： 原因背景影响解决方法测序化学反应的效率降低测序过程依赖一系列化学反应，例如在DNA测序中，DNA聚合酶需要将带有标记的核苷酸逐个添加到新生链上。随着时间的推移，反应体系中的化学物质浓度会逐渐降低，导致反应效率下降。这使得一些碱基在延伸过程中未能正确地被识别或添加，从而增加错误的概率。在测序过程中，使用高浓度的反应试剂和优化的反应条件，可以一定程度上减缓反应效率降低的速度。信号检测的难度增加在测序过程中，每个碱基的识别依赖于信号的检测。对于长读长的测序，信号在测序过程的后期可能会变得更加微弱或模糊。信号的微弱和模糊使得准确判断碱基类型变得更加困难，从而导致质量值下降。例如，在荧光标记测序中，荧光信号可能会因为光漂白或背景噪声而逐渐减弱。采用先进的光学系统和信号处理算法，可以提高信号检测的灵敏度和准确性，进而提高长读长测序的质量值。背景噪声的积累在测序过程中，背景噪声不可避免地存在。这些噪声可能是由于仪器的电子噪声、化学反应中的非特异性信号等引起的。随着测序读长的增加，背景噪声的影响会逐渐积累，使得信噪比降低。这使得在后期的测序过程中，正确识别碱基变得更加困难。通过优化测序平台的硬件设计，提高仪器的稳定性和信号处理能力，可以降低背景噪声的影响。测序仪的校准问题测序仪在开始测序时通常会进行校准，以确保信号检测的准确性。随着测序过程的进行，仪器的性能可能会发生微小的变化，导致校准参数不再完全准确。这会影响后期碱基的识别和质量值。定期对测序仪进行维护和校准，可以减缓这种影响。文库制备和模板质量测序的质量也受到文库制备和模板质量的影响。文库制备过程中存在错误或模板不纯，可能会导致测序过程中的错误积累。在长读长测序中，这些错误可能会更加明显。通过优化文库制备过程，提高模板的纯度和质量，可以提高测序的整体质量。 1.3 质量热图 当使用Illumina库，并且数据保留了原来的序列标识符，可输出下图。其中编码的是flowcell tile，可以知道rea的tile来源以及每个tile的质量分数，以分析质量差是否只与flowcell的某一些title区域相关。图中显示了每个title的平均质量的偏差。冷（blue）颜色质量值好（左图），热（red）颜色表明质量差（右图）。 1.4 密度曲线 每条序列的平均质量值的密度曲线图，好的数据质量呈现一个偏右拖尾分布，表示绝大分read平均质量值都落到Q值高的部分；差的数据质量呈现偏左拖尾及多峰分布。 1.5 ATGC碱基分布 ATGC碱基分布图：横轴表示各碱基位置，纵轴表示碱基百分比。碱基类型分布检查用于检测有无AT、GC分离现象，由于RNA-Seq所测的序列为随机打断的cDNA片段，因随机性打断及碱基互补配对原则，理论上，G和C、A和T的含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线。由于Reads 5’端的前几个碱基为随机引物序列，存在一定的偏好性，因此会在碱基分布图中出现前端波动较大的现象，属于正常情况。 1.6 GC含量分布 GC含量分布图：此模块测量文件中每个序列的整个长度的GC内容，并将其与建模的正态分布的GC内容进行比较。横轴为平均GC含量，纵轴为每个GC含量对应的序列数量；蓝线为理论分布，红线为测量值，二者越接近越好。 GC含量异常时，可以排查数据中是否存在大量的核糖体rna序列，核糖体rna序列过多意味着测到的数据中有效数据含量少，不能满足后续分析。 1.7 N含量分布 N含量分布图：N表示不能识别为ATGC中的任何一个碱基，为未知状态，碱基质量值Q值为0，N含量越多，说明数据质量越差。 1.8 read长度分布 read长度分布图：表示测序read的长度分布密度曲线图，横坐标为read长度，纵坐标为对应长度的read数。此图体现read 长度的分布范围，rawdata 一般所有序列都是一个固定长度，过滤后的cleandata 的read 长度较短占比过高时应该考虑rawdata中数据是否存在大量接头，导致过滤后数据短片段占比过高，此时存在文库构建失败问题。 1.9 序列重复性分布 序列重复性分布图：横坐标是duplication的次数，纵坐标是duplicated reads的百分比。转录组数据本身存在一定的read重复，qc报告非常容易出现警告或者红色xx，应该看具体的重复率，但重复率不应该过高比如超过70%。重复率过高提示该样本文库多样性少，即测得的read种类少，获得的基因序列有效信息少，可能影响后续定量准确性。 1.10 异常过表达序列 过表达reads可能存在的情况：

• 原始数据存在接头的可能性高，可查看过滤后的数据qc报告判断
• 来自高表达基因的片段，具有生物学意义
• 来自污染序列，导致文库中的这个序列比例异常高
• 文库检测到的序列多样性不高，重复率高

序列异常排查方向：

• 随机抽取fq文件中2000条read，与NCBI数据库做blast比对，判断是否有污染其他物种
• 与过表达的read序列做blast比对，查看其来源
• 过表达的read序列与文库构建的接头对比，看是否是接头

1.11 MultiQc整合FastQC结果 当测序文件过多时，使用MultiQc整合FastQC结果，得到一个所有样本的qc质控 xxxx.html 文件报告，可下载到本地使用浏览器打开查看。 multiqc *.zip -o ./ -n qc_fastqc 二、数据过滤 测序得到的原始序列含有接头序列或低质量序列，为了保证信息分析的准确性，需要对原始数据进行质量控制，得到高质量序列（Clean Reads）。 质控标准：

• 去除含接头的reads；
• 过滤去除低质量值数据，确保数据质量；
• 去除含有N（无法确定碱基信息）的比例大于5%的reads。

1.1 Trim Galore [Trim Galore](Babraham Bioinformatics - Trim Galore![2]) 围绕 Cutadapt[3] 和 FastQC[4] 的包装脚本工具，可始终如一地将质量和适配器修剪应用于 FastQ 文件，并为 MspI 消化的 RRBS 型 （Reduced Representation Bisufite-Seq） 文库提供一些额外的功能。 参数描述-j/--cores使用线程数，默认为1-q/--quality切除质量得分低于设置值的序列，默认值20--phred33/--phred64使用不同质量得分作为Phred得分标准，默认选项是phred33-a/--adapter输入adapter序列，可选，软件会自动搜索可能性最高的平台对应的adapter。自动搜选的平台有三个，分别是：--illumina，--nextera，--small_rna--length长度小于设定值的reads将被丢弃，Default: 20 bp.--max_length长度大于设定值的reads将被丢弃--stringency限定最少与adaptor序列重叠的碱基数，默认接头序列：AGATCGGAAGAGC--paired对于双端测序的一对reads，如果其中一个不合格被剔除，另一个也将被剔除-o/--output_dir设定输出目录，必须是已存在的目录，否则运行将报错。--fastqc当剪切结束后用默认选项对结果文件进行fastqc分析--max_n去除含有N碱基数大于n的序列

使用trim_galore软件试运行过滤一个fq文件

# 激活小环境
conda activate rna

# 进入到 cleandata 目录
cd$HOME/project/Human-16-Asthma-Trans/data/cleandata/

# trim
# -j：使用cpu数3
# -q：过滤20以下低质量Q值
# --length 20：过滤最小read长度，read被剪切掉接头后长度会发生变化，小于这个长度的会被过滤
# --max_n 3：read中N含量最大值，大于这个数的read会被过滤，N默认为 read长度的5%，这里N=63*0.05~=3
# --fastqc：对过滤后的cleandata 再运行一次质量评估 fastqc 
# --stringency：read与接头重叠的碱基个数
# -o：输出目录
trim_galore -j 3 -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ./ ../rawdata/SRR1039510_1.fastq.gz ../rawdata/SRR1039510_2.fastq.gz

多线程后台运行

# 激活小环境
conda activate rna

# 进入到 cleandata 目录
cd$HOME/project/Human-16-Asthma-Trans/data/cleandata/

# 先生成一个变量,为样本ID
# awk -F'/' '{print $NF}' 指定'/'分隔符取最后字段
ls ../rawdata/*_1.fastq.gz | awk -F'/''{print $NF}' | cut -d'_' -f1 >ID

# 查看ID类容
cat ID
# ID内容如下：
# SRR1039510
# SRR1039511
# SRR1039512

# -j:线程数，这里使用6个cpu线程
cat ID | whileread id
do
echo"trim_galore -j 3 -q 20 --length 20 --max_n 3 --stringency 3 --fastqc --paired -o ./  ../rawdata/${id}_1.fastq.gz ../rawdata/${id}_2.fastq.gz"
done >trim.sh

提交 trim.sh 任务至后台，可一次往服务器上面投递多个任务。

技巧：

• 用ParaFly并行任务，给你的分析提提速！
• 玩转服务器
• 玩崩服务器

# 安装
conda activate rna
conda install parafly -y

# 运行
# -CPU：这里一次投递3个样本
# 1>trim.log  2>&1 生成任务日志
# nohup …… & 提交后台命令
nohup ParaFly -c trim.sh -CPU 3 1>trim.log  2>&1 &
# 需要注意使用的总cpu数量，total=单个任务运行的线程数*一次并行的任务个数。
# 这里total = 3*3，即使用了9个线程。
# 使用htop看服务器的总线程数，最好一个任务不要超过总线程数的1/3-1/2。

# 使用MultiQc整合FastQC结果
multiqc *.zip -n qc_trim

常见报错：

• 库文件版本报错，导致命令调用失败

• 没有激活正确的环境，导致当前环境找不到命令

如何结束任务进程？

• 终止任务Ctrl+C只对直接运行的任务有效，后台任务无效
• kill命令kill -9 %numnum为jobs命令返回的ID；kill -9 PID采用ps fx或htop获取PID任务编号
• 暂停任务Ctrl+Z只对直接运行的任务有效，后台任务无效，可用fg %num将任务调至前台，用bg %num将任务调至后台。

1.2 参考基因组准备

1. 基因组文件：fasta
2. 注释文件：gff3/gtf

常用参考基因组数据库

• [√] Ensembl[5]
• [ ] NCBI[6]
• [ ] UCSC[7]

• Common name：物种的通俗名
• Scientific name：科学用名，一般使用斜体的拉丁字表示
• Taxon ID：物种ID号， NCBI 数据库中使用频繁

• ensembl 物种参考基因组版本号：物种为人的版本GRCh38，大的版本号，后面还有不同的小版本如release-112，不同的reasele。
• NCBI数据库的版本物种号：GCA_000001635.9
• UCSC：hg19/GRCh37, hg38/GRCh38, human genome

下载fasta、GTF、GFF3文件

参考基因组的版本不需要频繁更新，1-2年更新一次即可。

## 参考基因组准备:注意参考基因组版本信息
# 下载，Ensembl：.html
# /

# 进入到参考基因组目录
mkdir -p $HOME/database/GRCh38.113
cd$HOME/database/GRCh38.113

# 下载基因组序列axel  curl   wget
# wget在服务器下载的命令 -c断点续传，continue 
nohup wget -c .GRCh38.dna.primary_assembly.fa.gz >dna.log &

# 下载转录组序列
nohup wget -c .GRCh38.113.chr.gtf.gz >rna.log &

# 下载基因组注释文件
nohup wget -c .GRCh38.113.chr.gtf.gz >gtf.log &

nohup wget -c .GRCh38.113.chr.gff3.gz >gff.log&

# 上述文件下载完整后，再解压；否则文件不完整就解压会报错
# 再次强调，一定要在文件下载完后再进行解压！！！
nohup gunzip *gz >unzip.log &

Ensembl数据库基因命名方式

• ENS[species prefix][feature type prefix][a unique eleven digit number]
• Feature prefixes: E(exon), FM(Ensemblproteinfamily), G(gene), GT(gene tree), P(protein), R(regulatory feature), T(transcript)
• mouse gene：ENSMUSG###########

参考资料

[1]

FastQC: /

[2]

Trim Galore: /

[3]

Cutadapt: /

[4]

FastQC: /

[5]

Ensembl: .html?redirect=no

[6]

NCBI: .shtml

[7]

UCSC: /