测序数据从下载、质控、比对到计数的上游分析流程和问题记录。

Overview: RNA-Seq WorkFlow

部分参考： 详解RNAseq数据全流程分析，从fastq文件到富集分析一站式解决

质量控制（Quality Control，QC）：

  • 目的：评估测序数据的质量，检测潜在问题。
  • 步骤：使用工具（如 FastQC）分析测序数据，检查读长分布、质量分布、GC含量、过度重复、接头含量等。

数据预处理：

  • 目的：去除低质量 reads、去除接头、进行质量修剪。
  • 步骤：使用工具（如 fastp 或 Trimmomatic）对原始测序数据进行预处理，提高后续分析的准确性。

比对（Alignment）：

  • 目的：将测序 reads 映射到参考基因组上。
  • 步骤：使用比对工具（如 HISAT2、STAR、TopHat）进行比对。生成 SAM/BAM 文件，其中包含每个 read 的映射位置。

排序和索引（Index）：

  • 目的：提高数据的检索速度。
  • 步骤：使用工具（如 samtools）对比对后的 BAM 文件进行排序，并为 BAM 文件创建索引index。

计数（Count）：

  • 目的：统计每个基因的 reads 数量，用于后续表达量分析。
  • 步骤：使用工具（如 featureCounts 或 HTSeq）进行基因计数。生成一个表达矩阵，描述每个样本中每个基因的 reads 数量。

表达量分析：

  • 目的：研究基因的表达水平，找到差异表达基因。
  • 步骤：使用差异表达分析工具（如 DESeq2、edgeR、limma）进行统计分析，识别在不同条件下表达水平显著差异的基因。

功能注释：

  • 目的：对差异表达基因进行生物学功能的注释。
  • 步骤：做（KEGG、GO、GSEA、GSEVA）功能富集分析，了解差异表达基因的生物学过程、分子功能和细胞组分。

可视化和解释：

  • 目的：将分析结果可视化，提高对实验结果的理解。
  • 步骤：使用图形工具（如 R 中的 ggplot2、pheatmap等）绘制差异表达基因的热图、散点图、火山图等，帮助解释实验结果。

排序是什么？在哪里排？给谁排？比对的时候不好像也排序了吗？ 可以看最后补充内容节What is index

Preparation: Tools & Environment Variables

Tools

  • SRA-tools安装：SRA-tools，是一个NCBI数据下载工具。
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
# 创建名为rs0112的RNA-seq虚拟环境  
conda create --name rs0112  
# 激活虚拟环境  
conda activate rs0112
# 安装不同步骤需要的软件：
# 质控  
conda install fastqc  
# 预处理  
conda install trimmomatic  
conda install fastp  
# 比对  
conda install hisat2  
# 将sam文件转成bam文件  
conda install samtools  
# 计算count值  
conda install subread
  • 在samtools的安装上，专题：RNA-Seq ww上游分析学习 - yuanj's blog说这里有坑，建议安装1.9版本的，但我没管，事实上，这在后面确实遇到了坑。但是采用了其他办法解决了。
  • subread的安装是为了使用其包含的featureCounts来计算count值。
  • 其他如有遗漏，随时安装到本项目环境即可

Environment variables

设置环境变量：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 把下面的写进.bashrc里便于调用
export RNA_HOME=/mnt/disk1/gmq/240113UpStream
export RNA_SRA=$RNA_HOME/sra

export RNA_FQ=$RNA_HOME/fastq_gz
export RNA_FQRE=$RNA_HOME/fastqcReport

export RNA_TRIM=$RNA_HOME/trim
export RNA_TRIMRE=$RNA_HOME/trimReport

export RNA_ISRR=$RNA_FQ/SRR2303755
export RNA_OSRR=$RNA_TRIM/SRR2303755

export RNA_CL=$RNA_HOME/clean/SRR2303755
export RNA_CLRE=$RNA_HOME/cleanReport

export RNA_GENOME=$RNA_HOME/genome/genome
export RNA_SAM=$RNA_HOME/sam
export RNA_BAM=$RNA_HOME/bam

准备文件夹：

1
2
3
4
5
6
7
8
9
10
mkdir -p $RNA_SRA
mkdir -p $RNA_FQ
mkdir -p $RNA_FQRE
mkdir -p $RNA_TRIM
mkdir -p $RNA_TRIMRE
mkdir -p $RNA_CL
mkdir -p $RNA_CL_RE
mkdir -p $RNA_SAM
mkdir -p $RNA_SAM
mkdir -p $RNA_HOME/genome

如有需要，随时创建文件夹，并添加需要的环境变量即可。

nohup

如果希望命令在后台运行，同时又希望在终端的交互，那么可以使用 nohup 结合 &。

1
nohup bash -c 'YOUR_COMMANDS' &

怎么知道程序确实在工作呢？

  • 一个直观的方法是看看文件大小有没有变化

当然使用tmux也是很方便的

Get Data: SRA-tools

参考链接： https://github.com/ncbi/sra-tools/wiki/08.-prefetch-and-fasterq-dump 拓展：vdb工具的使用

prefetch

下载sra文件到SRA文件夹：

1
nohup prefetch --option-file -O $RNA_SRA $RNA_HOME/SRR_list.txt

这里的SRR_list.txt就是SRR那6个名字，可以直接手打。或者使用GNU的parallel工具

1
2
3
4
parallel -j 1 echo {1}_{2} ::: BORED EXCITED ::: 1 2 3 > ids
# 所以也可以用这样的：
parallel -j 1 echo {1}{2} ::: SRR2303755 ::: 1 2 3 4 5 6 > myids
# 但反正手打粘贴也是一种快，就是不那么优雅。

-j 1: 指定并行执行的最大任务数为1，表示每次只执行一个任务。

GNU 是什么？

  • GNU代表“GNU's Not Unix（GNU不是Unix）”。它是一个自由软件运动的产物，旨在创建一个类Unix操作系统的替代品，但是完全是自由开源的。由理查德·斯托曼（Richard Stallman）于1983年创始，并由自由软件基金会（Free Software Foundation）支持。
  • GNU计划的目标是提供一个完全自由、开源、类Unix的操作系统。虽然Linux内核与GNU工具一起被广泛用于创建现代GNU/Linux操作系统，但严格来说，GNU/Linux指的是使用Linux内核并包括GNU工具的操作系统。
  • GNU项目产生了许多重要的自由软件工具和组件，如GCC（GNU编译器集）、Emacs文本编辑器、Bash shell等。GNU的哲学强调用户自由，包括运行、修改、分享和贡献软件的自由。

fastq-dump/fasterq-dump

使用 fastq-dump 将其拆分：

1
fastq-dump --gzip --split-3 -O $RNA_FQ $RNA_SRA/SRR*.sra
  • –gzip 是将拆分的 fastq 文件压缩归档为 gz 格式
  • –split-3 是将文件拆分为正向序列和逆向序列 最后会得到1个.fastq.gz，1个正向_1.fastq.gz，1个逆向_2.fastq.gz 然后把他们挪到事先准备好的fastq_gz文件夹里面备用。

Prepare reference genome

准备参考基因组：

Reference Genomes | Griffith Lab (rnabio.org)

这么多数据库如何找参考基因组？应该以哪个为准？

  • 目前来看，接触到的教程都在使用Ensembel，将NCBI提供的ID转换为Ensembl ID。可以使用生物信息学工具如 biomaRt 包来进行这种映射。

我首先找到的是NCBI提供的Mus musculus genome assembly GRCm39 - NCBI - NLM (nih.gov)，因为给出的GEO数据库就是NCBI的，在给的GEO数据库页面就能找到：BioProject>>genome>>reference genome>>curl下载。下好的是zip文件，但是它并没有任何后缀，重命名加上.zip后缀，然后再解压缩，参考基因组数据存放在data文件夹内。

1
2
conda install unzip
unzip [YOUR_ZIP_FILE]

另外把这个geo的下载方法也备注于此：使用wget批量下载geo数据集的全部文件 | 生信菜鸟团 (bio-info-trainee.com)

1
wget -O filename.tar -o download.log -b URL

下载geo链接先找到文件夹所在位置，通过 download FAMILY --> Parent... --> supl --> copy 数据链接格式： https://ftp.ncbi.nlm.nih.gov/geo/series/GSE207nnn/GSE207502/suppl/GSE207502_RAW.tar

1
wget -O GSE244515_RAW.tar -o download.log -b URL https://ftp.ncbi.nlm.nih.gov/geo/series/GSE244nnn/GSE244515/suppl/GSE244515_RAW.tar

做到最后查看作者给出的DE结果，使用了Ensembl ID，因此，现在面临两条路：1. 重新用Ensembl数据库的参考基因组再次从alignment开始重新进行上游处理，2. 使用biomaRt进行ID映射。

首先解决第一个问题：如何下载Ensembl的参考基因组？ FTP Download (ensembl.org) 这里提到：Custom data sets: If you want to filter or customise your download, please try Biomart, a web-based querying tool.先不管这个biomaRt工具，先使用wget工具下载：

The complete data from which these files were obtained can be found at: ftp://ftp.ensembl.org/pub/release-86/fasta/homo_sapiens/dna/. You could use wget to download the Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz file, then unzip/untar. Reference Genomes | Griffith Lab (rnabio.org)

所以，我使用以下方法下载。

1
2
3
4
5
6
cd $RNA_HOME
mkdir ensembl_data
cd ensembl_data

# 下载参考基因组.fa
$ nohup wget https://ftp.ensembl.org/pub/current_fasta/mus_musculus/dna/Mus_musculus.GRCm39.dna.toplevel.fa.gz &

ENSEMBL FTP SITE: Based on Ensembl annotations only. Available for many species. http://useast.ensembl.org/info/data/ftp/index.html

1
2
# 下载注释.gtf
$ nohup wget https://ftp.ensembl.org/pub/current_gtf/mus_musculus/Mus_musculus.GRCm39.111.gtf.gz &

当然，这里下载的时候可以下载gff3，gtf等，我这里下载了gft。

使用biomaRt进行ID映射

Quality Control: FastQC

1
2
cd $RNA_FQ # 注意当前工作路径
nohup bash -c 'fastqc SRR*.fastq.gz -o $RNA_FQRE' &

报告的解读为下一步进行数据预处理提供指导： FastQC - Research Computing Documentation (usf.edu) FastQC aims to provide a QC report which can spot problems which originate either in the sequencer or in the starting library material. FastQC 的目的是提供一份质量控制报告，发现源于测序仪或起始文库材料的问题。

FastQC can be run in one of two modes. It can either run as a stand alone interactive application for the immediate analysis of small numbers of FastQ files, or it can be run in a non-interactive mode where it would be suitable for integrating into a larger analysis pipeline for the systematic processing of large numbers of files. FastQC 可以两种模式之一运行。它既可以作为一个独立的交互式应用程序运行，对少量 FastQ 文件进行即时分析，也可以在非交互模式下运行，适合集成到更大的分析管道中，对大量文件进行系统处理。

下面是选几个进行翻译总结，其他的用到再看吧

Per Base Sequence Quality

每碱基测序质量 Per Base Sequence Quality (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 长时间运行后质量普遍下降以及接头读取问题：根据平均质量和接头来综合地进行质量修剪
  2. 在运行的早期出现了短暂的质量损失：查看每层质量图(per tile quality plot)发现，不能做质量修剪，但需要在后续的比对和组装过程考虑屏蔽碱基(masking bases)
  3. 读长不同：查看读长分布模块，确定有多少序列触发了错误

Common reasons for warnings

The most common reason for warnings and failures in this module is a general degradation of quality over the duration of long runs. In general sequencing chemistry degrades with increasing read length and for long runs you may find that the general quality of the run falls to a level where a warning or error is triggered. 本模块中出现警告和故障的最常见原因是长时间运行后质量普遍下降。一般来说，测序化学性质会随着读取长度的增加而下降，对于长时间运行，您可能会发现运行的总体质量下降到触发警告或错误的水平。

If the quality of the library falls to a low level then the most common remedy is to perform quality trimming where reads are truncated based on their average quality. 如果文库质量较低，最常见的补救措施就是进行质量修剪，根据平均质量截断读数。

For most libraries where this type of degradation has occurred you will often be simultaneously running into the issue of adapter read-through so a combined adapter and quality trimming step is often employed. 对于大多数发生了这种类型降解的文库来说，你往往会同时遇到接头读穿的问题，因此通常会采用一个综合的接头和质量修剪步骤。

接头读取（adapter read-through）是什么？ 就是接头(adapter)序列也被读了的意思，所以要去接头

Another possibility is that a warn / error is triggered because of a short loss of quality earlier in the run, which then recovers to produce later good quality sequence. 另一种可能的情况是，由于在运行的早期出现了短暂的质量损失，从而触发了警告/错误，但随后又恢复并产生了质量良好的序列。

This can happen if there is a transient problem with the run (bubbles passing through a flowcell for example). You can normally see this type of error by looking at the per-tile quality plot (if available for your platform). 如果运行过程中出现瞬时问题（例如气泡通过流动池），就会出现这种情况。通常情况下，您可以通过查看每层质量图（如果您的平台有的话）来发现这种类型的错误。

每层质量图(per-tile quality plot)是什么？ 也是报告中的一部分，就是那个全是蓝色的那张图，蓝色代表很好的意思，绿色代表中等质量，红色代表低质量。Per Tile Sequence Quality (babraham.ac.uk)这个

In these cases trimming is not advisable as it will remove later good sequence, but you might want to consider masking bases during subsequent mapping or assembly. 在这种情况下，修剪是不可取的，因为它会删除后面的良好序列，但您可能需要考虑在后续的映射或组装过程中屏蔽碱基(masking bases)

屏蔽碱基(masking bases)是什么？ 在测序数据处理的早期阶段，可以使用工具（例如FastQC）进行质量控制和过滤。这些工具通常会标记低质量的碱基，你可以选择在后续的分析中排除或掩码这些碱基。根据per-base sequence quality图表中的质量分数信息，你可以确定要掩码的低质量碱基的位置。一般来说，质量分数较低的碱基可能需要进行掩码。

  1. 使用其他工具（例如sed、awk等）或脚本语言（例如Python或Perl）来在原始序列文件中标记或替换低质量碱基。这可能涉及到将低质量碱基替换为特殊字符，或者在需要掩码的位置添加标记。
  2. 请注意，这种方法是一种手动的方式，具体的步骤和操作会根据你的数据和需求而有所不同。在进行任何碱基掩码的操作之前，请确保了解你的数据和分析目的，以避免潜在的错误或误导性的结果。

If your library has reads of varying length then you can find a warning or error is triggered from this module because of very low coverage for a given base range. Before committing to any action, check how many sequences were responsible for triggering an error by looking at the sequence length distribution module results. 如果您的文库有不同长度的读数，那么您会发现该模块会触发警告或错误，因为特定碱基范围的覆盖率非常低。在采取任何行动之前，请查看序列长度分布模块的结果，检查有多少序列触发了错误。

Per Tile Sequence Quality

Summary

This graph will only appear in your analysis results if you're using an Illumina library which retains its original sequence identifiers. Encoded in these is the flowcell tile from which each read came. The graph allows you to look at the quality scores from each tile across all of your bases to see if there was a loss in quality associated with only one part of the flowcell. 只有在使用保留了原始序列标识符的 Illumina 文库时，该图才会出现在分析结果中。其中编码的是每个读数所来自的流式细胞片段。通过该图，您可以查看来自所有碱基的每个片段的质量得分，以了解是否只有流式细胞的一部分出现了质量损失。

The plot shows the deviation from the average quality for each tile. The colours are on a cold to hot scale, with cold colours being positions where the quality was at or above the average for that base in the run, and hotter colours indicate that a tile had worse qualities than other tiles for that base. In the example below you can see that certain tiles show consistently poor quality. A good plot should be blue all over. 该图显示了每个tile的质量与平均质量的偏差。颜色由冷到热，冷色代表在运行中质量处于或高于该基数平均水平的位置，而较热的颜色表示的质量比该基数的其他tile差。在下面的示例中，您可以看到某些tile的质量一直很差。好的地块应该全部为蓝色。

Reasons for seeing warnings or errors on this plot could be transient problems such as bubbles going through the flowcell, or they could be more permanent problems such as smudges on the flowcell or debris inside the flowcell lane. 在此图上看到警告或错误的原因可能是瞬时问题，如气泡穿过流量计，也可能是更持久的问题，如流量计上的污点或流量计通道内的碎屑。

Common reasons for warnings

Whilst warnings in this module can be triggered by individual specific events we have also observed that greater variation in the phred scores attributed to tiles can also appear when a flowcell is generally overloaded. In this case events appear all over the flowcell rather than being confined to a specific area or range of cycles. We would generally ignore errors which mildly affected a small number of tiles for only 1 or 2 cycles, but would pursue larger effects which showed high deviation in scores, or which persisted for several cycles. 虽然本模块中的警告可以由个别特定事件触发，但我们也观察到，当流室普遍超载时，瓦片的 phred 分数也会出现较大差异。在这种情况下，事件会出现在整个流动池中，而不是局限于某个特定区域或周期范围。我们通常会忽略仅在 1 或 2 个周期内对少量tile造成轻微影响的错误，但会关注分数偏差较大或持续数个周期的较大影响。

Per Sequence Quality Scores

Per Sequence Quality Scores (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 此处的警告通常表示运行中质量的普遍下降
  2. 可以通过修剪缓解此问题
  3. 若出现双模态(bi-modal)或复杂分布(complex distribution)，则结合每层质量进行评估

Common reasons for warnings

This module is generally fairly robust and errors here usually indicate a general loss of quality within a run. For long runs this may be alleviated through quality trimming. If a bi-modal, or complex distribution is seen then the results should be evaluated in concert with the per-tile qualities (if available) since this might indicate the reason for the loss in quality of a subset of sequences. 该模块通常相当稳健，此处的错误通常表示运行中质量的普遍下降。对于长时间运行，可以通过质量修剪来缓解这种情况。如果出现双模态或复杂分布，则应结合每层质量（如果有的话）对结果进行评估，因为这可能表明序列子集质量下降的原因。

Per Base Sequence Content

Per Base Sequence Content (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 样本中存在adapt dimmer or rRNA：样本问题，不知道怎么处理
  2. Biased Fragmentation：前12个碱基有问题，不需要修剪，不影响结果
  3. Biased composition library：特殊处理的文库，是正常的，不知道要不要处理
  4. aggressivley adapter trimmed：末段有偏差，假的不用管

Summary

Per Base Sequence Content plots out the proportion of each base position in a file for which each of the four normal DNA bases has been called. 每个碱基序列含量绘制出文件中四个正常 DNA 碱基中每个碱基位置被调用的比例。

It's worth noting that some types of library will always produce biased sequence composition, normally at the start of the read. 值得注意的是，某些类型的文库总是会产生有偏差的序列组成，通常是在读取开始时。

Libraries produced by priming using random hexamers (including nearly all RNA-Seq libraries) and those which were fragmented using transposases inherit an intrinsic bias in the positions at which reads start. 使用随机六聚体打底生成的文库（包括几乎所有的 RNA-Seq 文库）和使用转座酶片段化的文库在读数起始位置上存在固有偏差。

This bias does not concern an absolute sequence, but instead provides enrichement of a number of different K-mers at the 5' end of the reads. Whilst this is a true technical bias, it isn't something which can be corrected by trimming and in most cases doesn't seem to adversely affect the downstream analysis. It will however produce a warning or error in this module. 这种偏差并不涉及绝对序列，而是在读数的 5' 端富集了一些不同的K-mers。虽然这是一种真正的技术偏差，但并不能通过修剪来纠正，而且在大多数情况下似乎不会对下游分析产生不利影响。不过，它会在本模块中产生警告或错误。

什么是Kmers？ K-mer是DNA、RNA或蛋白质序列中的连续k个碱基或氨基酸的片段。具体来说，对于DNA或RNA序列，K-mer就是长度为k的连续碱基序列。例如，在DNA序列 "ATCGATCG" 中，当k=3时，可能的3-mer包括 "ATC"、"TCG"、"CGA" 等。

Common reasons for warnings

  1. Overrepresented sequences: If there is any evidence of overrepresented sequences such as adapter dimers or rRNA in a sample then these sequences may bias the overall composition and their sequence will emerge from this plot. 代表性过高的序列： 如果有任何证据表明样本中存在过度代表的序列（overrepresented sequences），如接头二聚体（adapter dimers）或 rRNA，那么这些序列可能会对总体组成产生偏差，它们的序列将从该图中显示出来。

什么是接头二聚体（adapter dimers）？ 由两个adapter连接形成的

  1. Biased fragmentation: Any library which is generated based on the ligation of random hexamers or through tagmentation should theoretically have good diversity through the sequence, but experience has shown that these libraries always have a selection bias in around the first 12bp of each run. This is due to a biased selection of random primers, but doesn't represent any individually biased sequences. Nearly all RNA-Seq libraries will fail this module because of this bias, but this is not a problem which can be fixed by processing, and it doesn't seem to adversely affect the ablity to measure expression.

片段选择偏差： 任何基于随机六聚体连接或通过标记产生的文库，理论上都应该具有良好的序列多样性，但经验表明，这些文库在每次运行的前 12bp 左右总是存在选择偏差。这是由于随机引物的选择偏差造成的，但并不代表任何有个别偏差的序列。几乎所有的 RNA-Seq 文库都会因为这种偏差而无法通过这一模块，但这不是一个可以通过处理来解决的问题，而且它似乎不会对表达测量能力产生不利影响。

  1. Biased composition libraries: Some libraries are inherently biased in their sequence composition. The most obvious example would be a library which has been treated with sodium bisulphite which will then have converted most of the cytosines to thymines, meaning that the base composition will be almost devoid of cytosines and will thus trigger an error, despite this being entirely normal for that type of library 有偏差的成分库： 有些文库的序列组成本身就有偏差。最明显的例子就是用亚硫酸氢钠处理过的文库，这种处理会将大部分胞嘧啶转化为胸腺嘧啶，这意味着碱基组成中几乎没有胞嘧啶，因此会引发错误，尽管这对这种类型的文库来说是完全正常的。

那这种文库序列本身就有偏差的会影响到基因表达的分析吗？如果会，应该如何解决？ 文库序列本身的偏差可能对基因表达分析产生影响，特别是在涉及到某些特殊文库制备方法时。以文库中的碱基组成偏差为例，例如通过亚硫酸钠处理的文库，其中大多数的胞嘧啶（cytosine）被转化为胸腺嘧啶（thymine）这种情况下，基因表达分析的结果可能受到一些偏差的影响。

  1. 适当的数据预处理： 在进行基因表达分析之前，进行适当的数据预处理步骤，例如质量控制、过滤低质量的读数、去除适配器等。这有助于减少一些制备偏差可能导致的分析误差。 2. 使用适合的工具和算法： 选择适合处理特殊文库偏差的分析工具和算法。有些工具可能对特殊情况具有更好的鲁棒性，能够更好地处理文库制备引入的偏差。 3. 进行偏差校正： 一些工具和软件提供了偏差校正的选项，可以根据已知的偏差信息进行校正。这可能包括在分析流程中纠正碱基组成的偏差。 4. 参考文献和文库制备方法的文档： 了解使用的文库制备方法并查阅相关文献，以更好地理解可能的偏差来源，并采取适当的步骤来解决或纠正这些偏差。
  1. If you are analysing a library which has been aggressivley adapter trimmed then you will naturally introduce a composition bias at the end of the reads as sequences which happen to match short stretches of adapter are removed, leaving only sequences which do not match. Sudden deviations in composition at the end of libraries which have undergone aggressive trimming are therefore likely to be spurious. 如果您分析的文库经过了积极的接头修剪（aggressivley adapter trimmed），那么您自然会在读数末端引入成分偏差，因为恰好与短段接头匹配的序列会被移除，只留下不匹配的序列。因此，经过积极修剪的文库末端突然出现的成分偏差很可能是假的。

Per Sequence GC Content

Per Sequence GC Content (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 异常分布可能表示偏差子集或污染
  2. 尖锐峰值：可能是特定污染物（例如接头二聚体）造成，可能会被高代表性序列模块发现
  3. 多峰值：可能是不同物种的污染

Summary

This module measures the GC content across the whole length of each sequence in a file and compares it to a modelled normal distribution of GC content. 该模块测量文件中每个序列全长的 GC 含量，并将其与 GC 含量的正态分布模型进行比较。

In a normal random library you would expect to see a roughly normal distribution of GC content where the central peak corresponds to the overall GC content of the underlying genome. Since we don't know the the GC content of the genome the modal GC content is calculated from the observed data and used to build a reference distribution. 在一个正常的随机文库中，你会发现 GC 含量大致呈正态分布，其中的中心峰值与底层基因组的总体 GC 含量相对应。由于我们不知道基因组的 GC 含量，因此会根据观察到的数据计算出 GC 含量的模态，并以此建立参考分布。

An unusually shaped distribution could indicate a contaminated library or some other kinds of biased subset. A normal distribution which is shifted indicates some systematic bias which is independent of base position. If there is a systematic bias which creates a shifted normal distribution then this won't be flagged as an error by the module since it doesn't know what your genome's GC content should be. 异常形状的分布可能表明文库受到了污染，或存在其他类型的偏差子集。偏移的正态分布表明存在与基位置无关的系统偏差。如果系统性偏差导致正态分布偏移，模块不会将其标记为错误，因为它不知道基因组的 GC 含量应该是多少。

Common reasons for warnings

Warnings in this module usually indicate a problem with the library. Sharp peaks on an otherwise smooth distribution are normally the result of a specific contaminant (adapter dimers for example), which may well be picked up by the overrepresented sequences module. Broader peaks may represent contamination with a different species. 该模块中的警告通常表明文库存在问题。在原本平滑的分布上出现尖锐的峰值，通常是特定污染物（例如接头二聚体）造成的，这很可能会被高代表性序列模块发现。多峰值可能代表不同物种的污染。

Duplicate Sequences

Duplicate Sequences (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 特定子集富集或低复杂度污染物的存在往往会在图中右侧产生尖峰，
  2. 如果峰值持续出现在蓝色迹线中，则表明存在大量不同的高重复序列，这可能表明存在污染集或非常严重的技术重复
  3. 技术重复（如PCR选择性扩增）与生物重复（随机结果）：无法区别二者，但可以通过随机码（random barcode）区分这两种重复。
  4. ==为识别低表达量转录本而有意过度测序，导致高表达转录本的大量重复。这时出现的警告是我们希望看到的结果==，
  5. 对于起始点受限文库，受限的起始点将会产生巨大的重复水平，这不应被视为一个问题，也不应通过重复数据删除来消除

Summary

In a diverse library most sequences will occur only once in the final set. A low level of duplication may indicate a very high level of coverage of the target sequence, but a high level of duplication is more likely to indicate some kind of enrichment bias (eg PCR over amplification). 在多样性文库中，大多数序列在最终集合中只会出现一次。低水平的重复可能表明目标序列的覆盖率很高，但高水平的重复更有可能表明某种富集偏差（如 PCR 过度扩增）。

This module counts the degree of duplication for every sequence in a library and creates a plot showing the relative number of sequences with different degrees of duplication. 该模块计算文库中每条序列的重复程度，并绘制图表显示不同重复程度序列的相对数量。

To cut down on the memory requirements for this module only sequences which first appear in the first 100,000 sequences in each file are analysed, but this should be enough to get a good impression for the duplication levels in the whole file. Each sequence is tracked to the end of the file to give a representative count of the overall duplication level. To cut down on the amount of information in the final plot any sequences with more than 10 duplicates are placed into grouped bins to give a clear impression of the overall duplication level without having to show each individual duplication value. 为减少该模块的内存需求，只分析每个文件中前 100,000 个序列中首次出现的序列，但这足以让人对整个文件的重复程度有一个良好的印象。每个序列都会被追踪到文件的末尾，以便对整体重复水平进行代表性统计。为了减少最终图谱中的信息量，任何重复次数超过 10 次的序列都会被放入分组箱中，以便在不显示单个重复值的情况下对整体重复程度有一个清晰的印象。

Because the duplication detection requires an exact sequence match over the whole length of the sequence, any reads over 75bp in length are truncated to 50bp for the purposes of this analysis. Even so, longer reads are more likely to contain sequencing errors which will artificially increase the observed diversity and will tend to underrepresent highly duplicated sequences. 由于重复检测需要序列全长的精确匹配，因此在本分析中，任何长度超过 75bp 的读数都会被截断为 50bp。即便如此，较长的读数也更有可能包含测序错误，这将人为地增加观察到的多样性，并倾向于低估高度重复序列的代表性。

The plot shows the proportion of the library which is made up of sequences in each of the different duplication level bins. There are two lines on the plot. The blue line takes the full sequence set and shows how its duplication levels are distributed. In the red plot the sequences are de-duplicated and the proportions shown are the proportions of the deduplicated set which come from different duplication levels in the original data. 图中显示的是不同重复水平区间的序列在文库中所占的比例。图中有两条线。蓝线是完整序列集，显示其重复水平的分布情况。在红色图中，序列被去重复，所显示的比例是去重复序列集中来自原始数据中不同重复水平的序列所占的比例。

In a properly diverse library most sequences should fall into the far left of the plot in both the red and blue lines. A general level of enrichment, indicating broad oversequencing in the library will tend to flatten the lines, lowering the low end and generally raising other categories. More specific enrichments of subsets, or the presence of low complexity contaminants will tend to produce spikes towards the right of the plot. 在一个具有适当多样性的文库中，大多数序列应属于图中红线和蓝线的最左端。一般的富集水平表明文库中存在广泛的过度序列，这往往会使线条变平，降低低端，并普遍提高其他类别。更具体的子集富集或低复杂度污染物的存在往往会在图中右侧产生尖峰。

什么是子类富集（specific enrichments of subsets）？ 在一个更大的文库中的某些特定子集合中，特定的序列或类别在数量上相对富集的现象。这种富集可能与实验设计、文库制备或其他实验条件相关。

These high duplication peaks will most often appear in the blue trace as they make up a high proportion of the original library, but usually disappear in the red trace as they make up an insignificant proportion of the deduplicated set. If peaks persist in the blue trace then this suggests that there are a large number of different highly duplicated sequences which might indicate either a contaminant set or a very severe technical duplication. 这些高重复峰通常会出现在蓝色迹线中，因为它们在原始文库中占很大比例，但通常会在红色迹线中消失，因为它们在重复序列集中所占比例很小。如果峰值持续出现在蓝色迹线中，则表明存在大量不同的高重复序列，这可能表明存在污染集或非常严重的技术重复。

The module also calculates an expected overall loss of sequence were the library to be deduplicated. This headline figure is shown at the top of the plot and gives a reasonable impression of the potential overall level of loss. 该模块还能计算出如果对文库进行重复数据删除，预计会造成的整体序列损失。这个标题数字显示在图表的顶部，为潜在的总体损失水平提供了一个合理的印象。

Common reasons for warnings

The underlying assumption of this module is of a diverse unenriched library. Any deviation from this assumption will naturally generate duplicates and can lead to warnings or errors from this module. 本模块的基本假设是建立一个多样化的、非富集的文库。任何偏离这一假设的行为都会自然产生重复，并可能导致本模块发出警告或出现错误。

In general there are two potential types of duplicate in a library, technical duplicates arising from PCR artefacts, or biological duplicates which are natural collisions where different copies of exactly the same sequence are randomly selected. From a sequence level there is no way to distinguish between these two types and both will be reported as duplicates here. 一般来说，文库中有两种潜在的重复类型，一种是由于 PCR 人为因素造成的技术重复，另一种是生物重复，即随机选择完全相同序列的不同拷贝的自然碰撞。从序列层面上看，这两种类型无法区分，因此在这里都将作为重复序列报告。

A warning or error in this module is simply a statement that you have exhausted the diversity in at least part of your library and are re-sequencing the same sequences. In a supposedly diverse library this would suggest that the diversity has been partially or completely exhausted and that you are therefore wasting sequencing capacity. However in some library types you will naturally tend to over-sequence parts of the library and therefore generate duplication and will therefore expect to see warnings or error from this module. 该模块中的警告或错误只是说明您已经详尽讨论了至少部分文库的多样性，并正在重新测序相同的序列。在所谓的多样性文库中，这表明多样性已经部分或完全详尽了，因此您正在浪费测序能力。不过，在某些类型的文库中，您会自然而然地倾向于对文库中的部分序列进行过度测序，从而产生重复，因此会希望从该模块中看到警告或错误。

在什么类型的文库中我们期望看到过度测序从而导致这一模块报错？ --后面这段解释了：为了观察低表达的转录本，进行了过度测序，因此导致了高表达的部分一起发生了过度测序，从而导致了较高的总体重复率。常在高重复率区间产生峰值。

In RNA-Seq libraries sequences from different transcripts will be present at wildly different levels in the starting population. In order to be able to observe lowly expressed transcripts it is therefore common to greatly over-sequence high expressed transcripts, and this will potentially create large set of duplicates. This will result in high overall duplication in this test, and will often produce peaks in the higher duplication bins. 在 RNA-Seq 文库中，不同转录本的序列在起始群体中的表达水平相差悬殊。因此，为了观察低表达的转录本，通常会对高表达的转录本进行大量的过度测序，这可能会产生大量的重复序列。这将导致该测试的总体重复率较高，并经常在较高的重复率区间产生峰值。

This duplication will come from physically connected regions, and an examination of the distribution of duplicates in a specific genomic region will allow the distinction between over-sequencing and general technical duplication, but these distinctions are not possible from raw fastq files. A similar situation can arise in highly enriched ChIP-Seq libraries although the duplication there is less pronounced. 这种重复将来自物理上相连的区域，对特定基因组区域重复分布的研究将有助于区分过度测序和一般技术性重复，但原始 fastq 文件无法进行这些区分。在高度富集的 ChIP-Seq 文库中也会出现类似的情况，不过这种情况下的重复没有那么明显。

Finally, if you have a library where the sequence start points are constrained (a library constructed around restriction sites for example, or an unfragmented small RNA library) then the constrained start sites will generate huge dupliction levels which should not be treated as a problem, nor removed by deduplication. In these types of library you should consider using a system such as random barcoding to allow the distinction of technical and biological duplicates. 最后，如果您有一个序列起始点受限的文库（如围绕限制位点构建的文库，或未破碎的小 RNA 文库），那么受限的起始点将会产生巨大的重复水平，这不应被视为一个问题，也不应通过重复数据删除来消除。在这类文库中，您应该考虑使用随机条形码等系统来区分技术重复和生物重复。

什么是围绕限制位点构建的文库，或未破碎的小 RNA 文库（a library constructed around restriction sites for example, or an unfragmented small RNA library）？

  1. 围绕限制位点构建的文库： 这是指在文库构建过程中，DNA或RNA样本是围绕特定的限制酶切位点进行构建的。限制酶是能够识别特定DNA序列并在其附近切割DNA链的酶。在构建限制位点文库时，通常会选择特定的限制酶切割DNA或RNA，形成特定长度的片段。这种文库的特点是片段的起始和终止点通常与限制位点相关。这样的文库构建方式可能导致在相同或相似位置存在大量的重复序列。
  2. 未破碎的小 RNA 文库： 这是指在构建小RNA文库时，不对小RNA进行碎裂或切割。相反，整个小RNA分子被直接用于文库构建。由于小RNA通常具有相对较短的长度，整个分子的使用可能导致在文库中存在大量的相同或高度相似的序列，尤其是在起始点。

使用随机条形码（random barcoding）等区分技术重复和生物重复是如何实现的？ 在高度复杂的文库中，重复的序列可能来自于技术性的重复，例如PCR扩增过程中的重复，也可能来自于生物学的重复，即相同的生物学实体在样本中的多次表示。为了区分这两种重复，可以使用随机条形码（random barcoding）的方法。 随机条形码是在文库构建过程中引入的一种标识序列，它们被随机添加到每个DNA片段或RNA分子上。每个样本中具有相同生物学起源的分子被分配相同的随机条形码。在后续的数据分析中，可以使用这些随机条形码来区分来自相同生物学实体的序列，从而准确计算生物学上的重复，而不是仅仅基于技术上的重复。

Overrepresented Sequences

Overrepresented Sequences (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 这个高代表序列具有很重要的生物学意义/文库污染
  2. 可能有存在末尾的高代表性序列被遗漏，这是一种Warning的假阴性。
  3. 读数截断导致的假阳性。

Summary

A normal high-throughput library will contain a diverse set of sequences, with no individual sequence making up a tiny fraction of the whole. Finding that a single sequence is very overrepresented in the set either means that it is highly biologically significant, or indicates that the library is contaminated, or not as diverse as you expected. 一个正常的高通量文库会包含多种多样的序列，没有一个序列只占整个文库的极小部分。如果发现单个序列在这组序列中的代表性非常高，要么说明该序列具有非常重要的生物学意义，要么说明文库受到了污染，或者并不像您预期的那样具有多样性。

This module lists all of the sequence which make up more than 0.1% of the total. To conserve memory only sequences which appear in the first 100,000 sequences are tracked to the end of the file. It is therefore possible that a sequence which is overrepresented but doesn't appear at the start of the file for some reason could be missed by this module. 该模块列出所有占总数 0.1%以上的序列。为了节省内存，只有出现在前 100,000 个序列中的序列才会被追踪到文件末尾。因此，有可能出现这样的情况，即某个序列的代表性过高，但由于某种原因没有出现在文件的开头，从而被该模块遗漏。

For each overrepresented sequence the program will look for matches in a database of common contaminants and will report the best hit it finds. Hits must be at least 20bp in length and have no more than 1 mismatch. Finding a hit doesn't necessarily mean that this is the source of the contamination, but may point you in the right direction. It's also worth pointing out that many adapter sequences are very similar to each other so you may get a hit reported which isn't technically correct, but which has very similar sequence to the actual match. 程序将在常见污染物数据库中查找与每个超标序列匹配的序列，并报告找到的最佳匹配序列。匹配序列的长度必须至少为 20bp，且错配不超过 1 个。找到一个匹配点并不一定意味着这就是污染源，但可以为您指明正确的方向。还值得指出的是，许多接头序列彼此非常相似，因此您可能会得到一个技术上并不正确，但与实际匹配序列非常相似的命中报告。

Because the duplication detection requires an exact sequence match over the whole length of the sequence any reads over 75bp in length are truncated to 50bp for the purposes of this analysis. Even so, longer reads are more likely to contain sequencing errors which will artificially increase the observed diversity and will tend to underrepresent highly duplicated sequences. 由于重复检测需要整个序列长度上的精确序列匹配，因此本分析将长度超过 75bp 的读数截断为 50bp。即便如此，较长的读数也更有可能包含测序错误，这将人为地增加观察到的多样性，并倾向于低估高度重复序列的代表性。

Common reasons for warnings

This module will often be triggered when used to analyse small RNA libraries where sequences are not subjected to random fragmentation, and the same sequence may natrually be present in a significant proportion of the library. 该模块通常会在分析小 RNA 文库时触发，因为小 RNA 文库中的序列不会被随机片段化，而且同一序列可能天然存在于相当大比例的文库中。

什么是小RNA文库？怎么判断我们的文库是不是小RNA文库？ 小RNA文库是一种用于分析小RNA的文库，其中小RNA通常指的是长度在20到200核苷酸之间的RNA分子。小RNA在细胞中起着调控基因表达、参与细胞信号传导等重要功能。构建小RNA文库的过程通常与其他RNA-Seq文库有所不同，因为小RNA通常不需要经过RNA的逆转录和放大步骤。 对于判断一个文库是否是小RNA文库，可以考虑以下几个方面：

  1. 文库构建方法： 小RNA文库通常使用专门设计的文库构建方法，其中包括适配器连接和cDNA合成等步骤。具体的文库构建方法可以查阅实验室或研究论文中的详细描述。
  2. 测序数据长度分布： 小RNA通常具有较短的长度，一般在20到200核苷酸之间。因此，通过观察测序数据的长度分布可以初步判断是否为小RNA文库。
  3. 实验目的： 根据实验的具体目的，例如是否关注小RNA的调控功能，也可以间接判断文库是否为小RNA文库。
  4. 如果在分析中触发了针对小RNA文库的特定模块，说明该分析工具对小RNA文库进行了一些特殊的处理或考虑，这也是一个指示文库类型的信号。

Adapter Content

Adapter Content (babraham.ac.uk)

Warning原因与解决方式概括：

  1. 可以发现序列末端建立的通读接头序列
  2. 但如果有高代表性序列（如接头二聚体）可能会看不出其他偏差
  3. 做接头修剪有助于解决此问题

Summary

The Kmer Content module will do a generic analysis of all of the Kmers in your library to find those which do not have even coverage through the length of your reads. This can find a number of different sources of bias in the library which can include the presence of read-through adapter sequences building up on the end of your sequences. Kmer Content 模块会对文库中的所有 Kmers 进行通用分析，找出那些在整个读数长度上覆盖率不均匀的 Kmers。这可以在文库中发现多种不同的偏差来源，其中包括在序列末端建立的通读接头序列。

You can however find that the presence of any overrepresented sequences in your library (such as adapter dimers) will cause the Kmer plot to be dominated by the Kmers these sequences contain, and that it's not always easy to see if there are other biases present in which you might be interested. 不过，您可能会发现，如果您的文库中存在任何代表性过高的序列（如接头二聚体），就会导致 Kmer 图被这些序列所包含的 Kmers 所主导，而且并不总是很容易看出是否存在您可能感兴趣的其他偏差。

One obvious class of sequences which you might want to analyse are adapter sequences. It is useful to know if your library contains a significant amount of adapter in order to be able to assess whether you need to adapter trim or not. Although the Kmer analysis can theoretically spot this kind of contamination it isn't always clear. This module therefore does a specific search for a set of separately defined Kmers and will give you a view of the total proportion of your library which contain these Kmers. A results trace will always be generated for all of the sequences present in the adapter config file so you can see the adapter content of your library, even if it's low. 您可能要分析的一类明显的序列是接头序列。了解您的文库中是否含有大量的接头，对于评估是否需要进行接头修剪非常有用。虽然 Kmer 分析理论上可以发现这种污染，但并不总是很清楚。因此，该模块会对一组单独定义的 Kmers 进行特定搜索，并为您提供包含这些 Kmers 的文库总比例。接头配置文件中的所有序列都会生成结果跟踪，因此即使接头含量很低，您也能看到文库中的接头含量。

The plot itself shows a cumulative percentage count of the proportion of your library which has seen each of the adapter sequences at each position. Once a sequence has been seen in a read it is counted as being present right through to the end of the read so the percentages you see will only increase as the read length goes on. 图谱本身显示的是您的文库中在每个位置出现过每个适配序列的比例的累积百分比计数。一旦一个序列在读数中出现，它就会一直被计算到读数的末尾，因此您看到的百分比只会随着读数长度的增加而增加。

Common reasons for warnings

Any library where a reasonable proportion of the insert sizes are shorter than the read length will trigger this module. This doesn't indicate a problem as such - just that the sequences will need to be adapter trimmed before proceeding with any downstream analysis. 任何插入长度比读取长度短的文库都会触发该模块。这并不表示有问题，只是在进行任何下游分析之前，需要对序列进行接头修剪。

不同的报告略有差异，那对于每一个样本进行数据预处理的标准应当是统一的还是个性的？ 对于每一个样本进行数据预处理的标准是一个需要慎重考虑的问题。在大多数情况下，虽然样本可能有一些差异，但在数据预处理的过程中仍然可以采用统一的标准，尤其是在以下情况下：

  1. 实验设计的一致性： 如果样本来自于相似的生物学条件，实验设计和处理过程都是一致的，那么采用统一的数据预处理标准通常是合理的。
  2. 技术平台的一致性： 如果使用相同的测序平台、实验方案和文库构建方法，样本之间的技术差异可能较小，因此可以考虑使用相同的数据预处理流程。
  3. 数据质量的一致性： 如果样本的数据质量特征相似，例如在测序深度、测序错误率、适配器污染等方面表现一致，那么可以采用相同的预处理标准。
  4. 分析的一致性需求： 如果最终的数据分析流程对于输入数据的一致性要求较高，采用统一的预处理标准可以降低分析中的不确定性。 然而，在一些情况下，可能需要根据样本的特定情况调整数据预处理的标准，例如：
  5. 样本特异性问题： 如果某个样本存在特定的问题，比如质量较差、适配器残留严重，可能需要针对该样本采用个性化的预处理策略。
  6. 实验设计的差异： 如果实验设计中存在一些样本之间的差异，比如处理条件的不同，可能需要根据这些差异调整数据预处理。
  7. 技术平台的差异： 如果不同样本使用了不同的测序平台或处理流程，可能需要针对性地调整预处理步骤。 在确定是否使用统一标准还是个性标准时，关键是要权衡样本之间的一致性和差异性，以及数据分析的要求。在任何情况下，始终建议对数据进行仔细的质量评估，并在预处理过程中记录和报告相关的处理步骤，以确保最终的数据质量和可靠性。

Preprocessing: trimmomatic & fastp

trimmomatic

whereis trimmomatic命令可以查看trimmomatic在系统中存放的位置，从而查阅adapter文件夹的位置，便于使用ILLUMINACLIP参数

参考代码：

  1. from: www.yuanj.top
1
trimmomatic PE -phred33 SRR23037551_1.fastq.gz SRR23037551_2.fastq.gz -summary Mus_musculus.summary -baseout SRR23037551.fastq.gz LEADING:3 TRAILING:3 SLIDINGWINDOW:5:20 HEADCROP:13 MINLEN:35
  1. from: https://bioinformatics.stackexchange.com/questions/10678/removing-nextera-transposase-adapters-cutadapt
1
java -jar /opt/software/Trimmomatic/0.39-Java-1.8/trimmomatic-0.39.jar PE -phred33 1004_R1.fastq.gz 1004_R2.fastq.gz ../QC_data/1004_R1.fastq.gz ../QC_data/1004.qcup_R1.fastq.gz ../QC_data/1004_R2.fastq.gz ../QC_data/1004.qcup_R2.fastq.gz ILLUMINACLIP:/opt/software/Trimmomatic/0.39-Java-1.8/adapters/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 MINLEN:36
  1. from: https://speciationgenomics.github.io/Trimmomatic/
1
java -jar ~/Biosoft/Trimmomatic038/Trimmomatic-0.38/trimmomatic-0.38.jar PE -phred33 test_7942raw_1.fq.gz test_7942raw_2.fq.gz ./trim_out/output_forward_paired.fq.gz ./trim_out/output_forward_unpaired.fq.gz ./trim_out/output_reverse_paired.fq.gz ./trim_out/output_reverse_unpaired.fq.gz ILLUMINACLIP:/home/raewyn/Biosoft/Trimmomatic038/Trimmomatic-0.38/adapters/TruSeq2-PE.fa:2:30:10 SLIDINGWINDOW:5:20LEADING:20 TRAILING:20 MINLEN:75

综抄后试运行代码（仅对第一条SRR）：

1
trimmomatic PE -phred33 $RNA_FQ/SRR23037551_1.fastq.gz $RNA_FQ/SRR23037551_2.fastq.gz $RNA_1P $RNA_1U $RNA_2P $RNA_2U ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35

运行：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
## COMMANDS
gmq@zhanglab /mnt/disk1/gmq/240113_UpStream/trim_out
$ trimmomatic PE -phred33 $RNA_FQ/SRR23037551_1.fastq.gz $RNA_FQ/SRR23037551_2.fastq.gz $RNA_1P $RNA_1U $RNA_2P $RNA_2U ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35

## RESULT
TrimmomaticPE: Started with arguments:
 -phred33 /mnt/disk1/gmq/240113_UpStream/fastq_gz/SRR23037551_1.fastq.gz /mnt/disk1/gmq/240113_UpStream/fastq_gz/SRR23037551_2.fastq.gz /mnt/disk1/gmq/240113_UpStream/trim_out/1p.fq.gz /mnt/disk1/gmq/240113_UpStream/trim_out/1u.fq.gz /mnt/disk1/gmq/240113_UpStream/trim_out/2p.fq.gz /mnt/disk1/gmq/240113_UpStream/trim_out/2u.fq.gz ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35
Using PrefixPair: 'AGATGTGTATAAGAGACAG' and 'AGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTCCGAGCCCACGAGAC'
Using Long Clipping Sequence: 'CTGTCTCTTATACACATCTGACGCTGCCGACGA'
ILLUMINACLIP: Using 1 prefix pairs, 4 forward/reverse sequences, 0 forward only sequences, 0 reverse only sequences
Input Read Pairs: 32696486 Both Surviving: 28380533 (86.80%) Forward Only Surviving: 2556324 (7.82%) Reverse Only Surviving: 684775 (2.09%) Dropped: 1074854 (3.29%)
TrimmomaticPE: Completed successfully

那么再看看这个运行之后的fastq_report 将trimmomatic处理的第一个文件(forward_paired)用fastqc处理，生成一个html文件和一个zip文件

1
fastqc $RNA_1P -O $RNA_HOME/trim_out_report

看起来么有问题，那么接下来就对所有的样本文件进行处理。

使用for循环来运行所有待处理的fastq文件有两种：

1
2
3
for ((i=1; i<=6; i++)); do 
  trimmomatic PE -phred33 "${RNA_ISRR}${i}_1.fastq.gz" "${RNA_ISRR}${i}_2.fastq.gz" -summary "${RNA_OSRR}${i}.summary" -baseout "${RNA_OSRR}${i}.fastq.gz" ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35 
done

或者：

1
for ((i=1; i<=6; i++)); do trimmomatic PE -phred33 "${RNA_ISRR}${i}_1.fastq.gz" "${RNA_ISRR}${i}_2.fastq.gz" -summary "${RNA_OSRR}${i}.summary" -baseout "${RNA_OSRR}${i}.fastq.gz" ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35; done

也就是说，在Bash脚本中，do，trimmomatic命令和done之间可以使用换行或分号进行分隔，但不能既不换行也不用分号，这样是无法运行的，会有>提示你输入未完成。而我们又要把所有命令用单引号包括在内，因此用分号写法就好了。

运行命令：

1
nohup bash -c 'for ((i=1; i<=6; i++)); do trimmomatic PE -phred33 "${RNA_ISRR}${i}_1.fastq.gz" "${RNA_ISRR}${i}_2.fastq.gz" -summary "${RNA_OSRR}${i}.summary" -baseout "${RNA_OSRR}${i}.fastq.gz" ILLUMINACLIP:/home/gmq/micromamba/envs/rs0112/share/trimmomatic-0.39-2/adapters/NexteraPE-PE.fa:2:30:10 SLIDINGWINDOW:5:20 LEADING:20 TRAILING:20 HEADCROP:13 MINLEN:35; done' &

然后对正序paired文件生成报告：

1
nohup bash -c 'for ((i=3; i<=6; i++)); do fastqc "${RNA_OSRR}${i}_1P.fastq.gz" -O "${RNA_TRIMRE}"; done' &

注意：一定要在文件名处加""，双引号中包裹了整个赋值语句，这样可以确保整个字符串被认为是一个变量，不然就会由于赋值错误而找不到文件。但是通配符 * 的行为与引号的使用是有区别的。当通配符在引号内时不会被展开，而认为*是文件名称的一部分

fastp

是一个更快的工具，用C++编写

下面给出一个搞笑bug：

1
2
# 这个是运行不了的
nohup bash -c 'for ((i=1; i<=6; i++)); do fastp -i "${RNA_ISRR}${i}_1.fastq.gz" -I "${RNA_ISRR}${i}_2.fastq.gz" -o "${RNA_CL}${i}_R1.fastq" -O "${RNA_CL}${i}_R2.fastq" -h "fastp_report${i}.html"; done' &
1
2
# 请运行这个
nohup bash -c 'for ((i=1; i<=6; i++)); do fastp -i "${RNA_ISRR}${i}_1.fastq.gz" -I "${RNA_ISRR}${i}_2.fastq.gz" -o "${RNA_CL}${i}_R1.fastq" -O "${RNA_CL}${i}_R2.fastq" -h "fastp_report${i}.html"; done' &

为什么呢？ 看起来报错的原因是由于命令中的空格符存在问题。在你的 bash -c 命令中，fastp 命令的参数之间有不正常的空格。在你的命令中，你使用了全角空格（fastp -i）而不是普通的空格。请确保在你的命令中使用标准的ASCII空格。

所以请小心全角空格，不然就会卡很久也不知道哪里有问题。

fastp的参考论文链接在：fastp: an ultra-fast all-in-one FASTQ preprocessor 这里面的背景部分提到：

A typical combination was the use of FASTQC for quality control, Cutadapt for adapter trimming and Trimmomatic for read pruning and filtering......In this paper, we present fastp, an ultra-fast tool to perform quality control, read filtering and base correction for FASTQ data. It includes most features of FASTQC + Cutadapt + Trimmomatic + AfterQC while running 2–5 times faster than any of them alone.

fastp和trimmomatic的比较（GPT）： fastp 和 Trimmomatic 都是用于对高通量测序数据进行质量控制和适配器修剪的工具，但它们之间有一些区别。以下是它们的一些主要区别：

  1. 速度：
  • fastp 是一个被设计为非常快速的工具，它通过并行处理来提高处理速度。它的设计目标之一是在保持质量的同时提供高性能。
  • Trimmomatic 也是一个流行的工具，但相对而言可能稍慢一些，特别是在处理大型数据集时。它使用的是Java编写，而且相对较旧，可能在某些情况下可能不如 fastp 那样高效。
  1. 适配器修剪算法：
  • fastp 使用了一种基于最小K-mer哈希值的适配器检测算法。它能够快速准确地检测适配器序列，并进行修剪。
  • Trimmomatic 也有适配器修剪功能，但它使用的是更传统的Seed-and-Extend算法，这在一些情况下可能相对慢一些。
  1. 功能和设置：
  • fastp 提供了一些附加功能，例如质量过滤、多种格式的输出报告、复杂的质量过滤算法等。它还支持多种测序平台的数据。
  • Trimmomatic 是一个功能强大的工具，支持多种修剪和过滤方式，包括LEADING、TRAILING、SLIDINGWINDOW等。它提供了详细的文档，允许用户进行灵活的参数设置。
  1. 平台支持：
  • fastp 可以处理 Illumina、Ion Torrent 和 PacBio 等不同平台的数据。
  • Trimmomatic 主要设计用于处理 Illumina 测序数据，虽然也可以用于其他平台。
  1. 语言：
  • fastp 使用C++编写。
  • Trimmomatic 使用Java编写。

fastp和trimmomatic的选择（GPT）： 以下是一些考虑因素：

  1. 速度： 如果你对处理速度有较高的要求，fastp 可能更适合，因为它被设计为非常快速的工具。
  2. 功能： 如果你需要更灵活的参数设置、更多的过滤和修剪选项，以及更复杂的报告输出，Trimmomatic 提供了更多的功能和选项。
  3. 适配器检测： fastp 使用了一种不同的适配器检测算法，可能在某些情况下更准确或更适用。
  4. 平台支持： 如果你处理多种平台的数据，fastp 提供了对不同平台的支持。
  5. 易用性： fastp 的命令行选项相对简单，可能更容易上手。Trimmomatic 的参数较多，可能需要更多的配置。

最终，选择使用哪个工具取决于你的具体需求和实验设计。

  • 你可以根据实验数据的性质、处理速度的要求、需要进行的过滤和修剪操作等因素来做出决策。
  • 最好的方法是在小规模数据上尝试这两个工具，看看它们在你的特定情境下的表现如何。

hisat

1
hisat2-build -p 4 /mnt/disk1/gmq/240113_UpStream/reference/genome.fna genome
1
2
3
4
5
6
7
8
9
10
11
# 解压下载好的基因组.fa.gz文件
nohup bash -c 'gzip Mus_musculus.GRCm39.dna.toplevel.fa.gz' &

# 解压下载好的注释.gtf.gz文件
nohup bash -c 'gunzip Mus_musculus.GRCm39.111.gtf.gz' &

# 使用hisat2构建index: engenome.1.ht2, engenome.2.ht2 ...
nohup bash -c 'hisat2-build -p 4 /mnt/disk1/gmq/240113UpStream/ensembl_data/engenome.fa engenome' &

# 把engenome基因组与index也移动到engenome文件夹中
nohup bash -c 'cp -t ../engenome engenome*' &
1
2
3
hisat2 -x genome/genome -1 trim_out/SRR2303755*_1P.fastq.gz -2 trim_out/SRR2303755*_2P.fastq.gz -S genome.sam

nohup bash -c 'hisat2 -x engenome/engenome -1 trim_out/SRR2303755*_1P.fastq.gz -2 trim_out/SRR2303755*_2P.fastq.gz -S engenome.sam' &
1
hisat2 -x genome/genome -1 trim_out/SRR2303755*_1P.fastq.gz -2 trim_out/SRR2303755*_2P.fastq.gz -S genome.sam
1
2
3
4
5
6
7
8
9
10
# 先复制过来再改
RNA_GENOME=$RNA_HOME/genome/genome.fna
/mnt/disk1/gmq/240113_UpStream/reference/
mkdir -p bam
parallel -j 1 echo {1}_{2} ::: BORED EXCITED ::: 1 2 3
cat ids | parallel "hisat2 -x $RNA_GENOME -1 trim_out/{}_1P.fastq.gz -2 trim_out/{}_2P.fastq.gz | samtools sort > bam/{}.bam"
# Index each BAM file.
cat ids | parallel "samtools index bam/{}.bam"
#  we can see the resulting files with:
ls -l bam

报错如果一直找不到genome文件，不要写.fna后缀

1
nohup bash -c 'hisat2 -x /mnt/disk1/gmq/240113UpStream/genome/genome -1 $RNA_TRIM/SRR23037551_1P.fastq.gz -2 $RNA_TRIM/SRR23037551_2P.fastq.gz -S $RNA_HOME/sam/SRR23037551.sam' &
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
## SRR23037551
bash: hisat2: command not found
28380533 reads; of these:
  28380533 (100.00%) were paired; of these:
    5948023 (20.96%) aligned concordantly 0 times
    20713019 (72.98%) aligned concordantly exactly 1 time
    1719491 (6.06%) aligned concordantly >1 times
    ----
    5948023 pairs aligned concordantly 0 times; of these:
      1721659 (28.95%) aligned discordantly 1 time
    ----
    4226364 pairs aligned 0 times concordantly or discordantly; of these:
      8452728 mates make up the pairs; of these:
        3792456 (44.87%) aligned 0 times
        2983353 (35.29%) aligned exactly 1 time
        1676919 (19.84%) aligned >1 times
93.32% overall alignment rate

用for循环：

1
2
3
4
5
6
7
8
9
10
11
nohup bash -c 'for ((i=2; i<=6; i++)) do hisat2 -x /mnt/disk1/gmq/240113UpStream/genome/genome -1 "${RNA_TRIM}/SRR2303755${i}_1P.fastq.gz" -2 "${RNA_TRIM}/SRR2303755${i}_2P.fastq.gz" -S "${RNA_HOME}/sam/SRR2303755${i}.sam"; done' &

nohup bash -c 'for ((i=1; i<=6; i++)) do hisat2 -x /mnt/disk1/gmq/240113UpStream/engenome/engenome -1 "${RNA_TRIM}/SRR2303755${i}_1P.fastq.gz" -2 "${RNA_TRIM}/SRR2303755${i}_2P.fastq.gz" -S "${RNA_HOME}/sam/ENSRR2303755${i}.sam"; done' &

## 在后面遇到了这样的报错：
[E::sam_parse1] SEQ and QUAL are of different length
[W::sam_read1] parse error at line 38488487
[bam_sort_core] truncated file. Aborting.
这可能是因为在运行的时候硬盘空间不足导致的，重新运行一次。

nohup bash -c 'hisat2 -x $RNA_HOME/engenome/engenome -1 "$RNA_TRIM/SRR23037556_1P.fastq.gz" -2 "$RNA_TRIM/SRR23037556_2P.fastq.gz" -S "$RNA_HOME/sam/ENSRR23037556.sam"' &

samtool

参考代码：专题：RNA-Seq 上游分析学习 - yuanj's blog

1
nohup bash -c 'samtools sort -n -@ 5 oryza_sativa.sam -o oryza_sativa' 
1
2
3
4
cd $RNA_BAM
nohup bash -c 'samtools sort -@ 4 $RNA_SAM/SRR23037551.sam -o $RNA_BAM/SRR23037551' &
## 报错如下：
samtools: error while loading shared libraries: libncurses.so.5: cannot open shared object file: No such file or directory

尝试安装libncurses解决问题：

1
2
3
4
5
6
7
8
9
$ mamba install -c conda-forge libncurses
## 报错如下
Pinned packages:
  - python 3.11.*

error    libmamba Could not solve for environment specs
    The following package could not be installed
    └─ libncurses does not exist (perhaps a typo or a missing channel).
critical libmamba Could not solve for environment specs

尝试安装libncurses=5解决问题：

1
2
3
4
5
6
7
8
9
$ mamba install -c conda-forge libncurses=5
## 报错如下：
Pinned packages:
  - python 3.11.*

error    libmamba Could not solve for environment specs
    The following package could not be installed
    └─ libncurses 5**  does not exist (perhaps a typo or a missing channel).
critical libmamba Could not solve for environment specs

重新开个环境rs0121

1
$ conda create -n rs0121 python=3.10

发现问题了，这里有个叫ncurses的软件，一查，发现它就是那个装不了的libncurses： ncurses is an open-source clone of the original Unix curses library. libcurses. usually points to libncurses.* to provide compatibility with the original library, so there would be no practical difference between using one over the other.* https://stackoverflow.com/questions/1517756/whats-the-difference-between-lcurses-and-lncurses-when-compiling-c-using-ncur

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
$ conda list -n rs0121
List of packages in environment: "/home/gmq/micromamba/envs/rs0121"

  Name              Version     Build               Channel
─────────────────────────────────────────────────────────────────
  _libgcc_mutex     0.1         conda_forge         conda-forge
  _openmp_mutex     4.5         2_gnu               conda-forge
  bzip2             1.0.8       hd590300_5          conda-forge
  ca-certificates   2023.11.17  hbcca054_0          conda-forge
  ld_impl_linux-64  2.40        h41732ed_0          conda-forge
  libffi            3.4.2       h7f98852_5          conda-forge
  libgcc-ng         13.2.0      h807b86a_3          conda-forge
  libgomp           13.2.0      h807b86a_3          conda-forge
  libnsl            2.0.1       hd590300_0          conda-forge
  libsqlite         3.44.2      h2797004_0          conda-forge
  libuuid           2.38.1      h0b41bf4_0          conda-forge
  libxcrypt         4.4.36      hd590300_1          conda-forge
  libzlib           1.2.13      hd590300_5          conda-forge
  ncurses           6.4         h59595ed_2          conda-forge
  openssl           3.2.0       hd590300_1          conda-forge
  pip               23.3.2      pyhd8ed1ab_0        conda-forge
  python            3.10.13     hd12c33a_1_cpython  conda-forge
  readline          8.2         h8228510_1          conda-forge
  setuptools        69.0.3      pyhd8ed1ab_0        conda-forge
  tk                8.6.13      noxft_h4845f30_101  conda-forge
  tzdata            2023d       h0c530f3_0          conda-forge
  wheel             0.42.0      pyhd8ed1ab_0        conda-forge
  xz                5.2.6       h166bdaf_0          conda-forge
(rs0112)

这重新看了一眼，确实是6.4版本没错了，那么就只需要把它降到5了

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
$ conda install ncurses=5
conda-forge/linux-64                                        Using cache
conda-forge/noarch                                          Using cache
bioconda/linux-64                                           Using cache
bioconda/noarch                                             Using cache

Pinned packages:
  - python 3.11.*


Transaction

  Prefix: /home/gmq/micromamba/envs/rs0112

  Updating specs:

   - ncurses=5


  Package        Version  Build               Channel           Size
──────────────────────────────────────────────────────────────────────
  Remove:
──────────────────────────────────────────────────────────────────────

  - hisat2         2.2.1  hdbdd923_6          bioconda        Cached
  - pip           23.3.2  pyhd8ed1ab_0        conda-forge     Cached
  - wheel         0.42.0  pyhd8ed1ab_0        conda-forge     Cached
  - setuptools    69.0.3  pyhd8ed1ab_0        conda-forge     Cached
  - python        3.11.0  h10a6764_1_cpython  conda-forge     Cached

  Downgrade:
──────────────────────────────────────────────────────────────────────

  - trimmomatic     0.39  hdfd78af_2          bioconda        Cached
  + trimmomatic     0.36  6                   bioconda         141kB
  - ncurses          6.4  h59595ed_2          conda-forge     Cached
  + ncurses          5.9  10                  conda-forge        1MB
  - readline         8.2  h8228510_1          conda-forge     Cached
  + readline         7.0  0                   conda-forge      859kB

  Summary:

  Remove: 5 packages
  Downgrade: 3 packages

  Total download: 2MB

──────────────────────────────────────────────────────────────────────


Confirm changes: [Y/n] y

这次再运行：

1
2
3
nohup bash -c 'samtools sort -@ 4 $RNA_SAM/SRR23037551.sam -o $RNA_BAM/SRR23037551' &
## nohup.out
[bam_sort_core] merging from 28 files...

featureCounts

1
nohup bash -c 'featureCounts -g ID -p -B -C -a $RNA_HOME/reference/features.gff -o counts.txt $RNA_BAM/SRR23037551.bam' &
1
2
3
ERROR: failed to find the gene identifier attribute in the 9th column of the provided GTF file.
The specified gene identifier attribute is 'gene_id' 
An example of attributes included in your GTF annotation is 'ID=exon-XR_004936710.1-1;Parent=rna-XR_004936710.1;Dbxref=GeneID:115487594,RFAM:RF00026,GenBank:XR_004936710.1,MGI:MGI:5455983;gbkey=ncRNA;gene=Gm26206;inference=COORDINATES: profile:INFERNAL:1.1.1;product=U6 spliceosomal RNA;transcript_id=XR_004936710.1'.

提示第9列不是gene_id？

1
2
3
4
5
6
7
8
9
10
11
$ head $RNA_HOME/reference/features.gff
##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
#!genome-build GRCm39
#!genome-build-accession NCBI_Assembly:GCF_000001635.27
#!annotation-date 04/05/2023
#!annotation-source NCBI RefSeq GCF_000001635.27-RS_2023_04
##sequence-region NC_000067.7 1 195154279
##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=10090
NC_000067.7     RefSeq  region  1       195154279       .       +       .      ID=NC_000067.7:1..195154279;Dbxref=taxon:10090;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA;strain=C57BL/6J

看到第9列的名字是ID

1
nohup bash -c 'featureCounts -g ID -p -B -C -a $RNA_HOME/reference/features.gff -o counts.txt $RNA_BAM/SRR23037551.bam' &
  • -p: 只能用在paired-end的情况中，会统计fragment而不统计read
  • -B: 在-p选择的条件下，只有两端read都比对上的fragment才会被统计
  • -C: 如果-C被设置，那融合的fragment（比对到不同染色体上的fragment）就不会被计数，这个只有在-p被设置的条件下使用
  • -T: 线程数目，1~32
  • -a: 基因组注释文件（gtf文件），也支持Gzipped文件格式
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
        ==========     _____ _    _ ____  _____  ______          _____  
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
	  v2.0.6

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                                                                            ||
||                           SRR23037551.bam                                  ||
||                                                                            ||
||             Output file : counts.txt                                       ||
||                 Summary : counts.txt.summary                               ||
||              Paired-end : yes                                              ||
||        Count read pairs : no                                               ||
||              Annotation : features.gff (GTF)                               ||
||      Dir for temp files : ./                                               ||
||                                                                            ||
||                 Threads : 1                                                ||
||                   Level : meta-feature level                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : not counted                                      ||
||   Min overlapping bases : 1                                                ||
||                                                                            ||
\\============================================================================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file features.gff ...                                      ||
||    Features : 1570915                                                      ||
||    Meta-features : 1570915                                                 ||
||    Chromosomes/contigs : 43                                                ||
||                                                                            ||
|| Process BAM file SRR23037551.bam...                                        ||
||    Paired-end reads are included.                                          ||
||    The reads are assigned on the single-end mode.                          ||
||    Total alignments : 67160967                                             ||
||    Successfully assigned alignments : 11721319 (17.5%)                     ||
||    Running time : 1.31 minutes                                             ||
||                                                                            ||
|| Write the final count table.                                               ||
|| Write the read assignment summary.                                         ||
||                                                                            ||
|| Summary of counting results can be found in file "counts.txt.summary"      ||
||                                                                            ||
\\============================================================================//
1
2
3
4
nohup bash -c 'featureCounts -g ID -p -B -C -a $RNA_HOME/reference/features.gff -o counts_2.txt $RNA_BAM/SRR23037552.bam' &

...

1
2
3
4
5
6
7
8
9
10
11
$ head $RNA_HOME/reference/counts_1.txt
# Program:featureCounts v2.0.6; Command:"featureCounts" "-g" "ID" "-p" "-B" "-C" "-a" "/mnt/disk1/gmq/240113UpStream/reference/features.gff" "-o" "counts.txt" "/mnt/disk1/gmq/240113UpStream/bam/SRR23037551.bam"
Geneid  Chr     Start   End     Strand  Length  /mnt/disk1/gmq/240113UpStream/bam/SRR23037551.bam
exon-XR_004936710.1-1   NC_000067.7     3172239 3172348 +       110     0
exon-XM_006495550.5-1   NC_000067.7     3740775 3741733 -       959     0
exon-XM_006495550.5-2   NC_000067.7     3491925 3492124 -       200     0
exon-XM_006495550.5-3   NC_000067.7     3283662 3287191 -       3530    0
exon-XM_006495550.5-4   NC_000067.7     3269956 3277540 -       7585    0
exon-NM_001011874.1-1   NC_000067.7     3740775 3741721 -       947     0
exon-NM_001011874.1-2   NC_000067.7     3491925 3492124 -       200     0
exon-NM_001011874.1-3   NC_000067.7     3284705 3287191 -       2487    0

Supplement

Background

CITE FROM：生物医学信息学_李兰娟_第八章_第二节_组学数据_转录组学数据分析

  1. 转录组数据的产生
  • 地位：RNA-seq技术是目前检测细胞和组织中全转录组数据的一种强有效的高通量测序方法。
  • 流程：测序的流程包括RNA样本制备、cDNA文库构建及高通量测序。
  • 作用：它能够在单核苷酸水平对任意物种的整体转录活动进行检测，在分析转录本的结构和表达水平的同时，还能发现未知转录本及未知基因组序列的物种，精确地识别可变剪切事件等，提供全面的转录组信息。
  • 优势：相对于传统的芯片杂交平台，转录组测序无需预先针对已知序列设计探针，并提供更精确的数字化信号，更高的检测通量以及更广泛的检测范围，是目前深入研究转录组复杂性的强大工具。
  • 应用：转录组数据可用于新基因和新转录本的预测，非编码RNA注释，基因的差异表达分析，SNP（单核苷酸多态）/INDEL（插入缺失）检测以及可变剪切事件检测等。
  1. 转录组数据分析 一个简单的流程：
  • 首先进行数据的比对，即将测序产生的reads进行质量控制，过滤掉低质量的reads后用比对软件将这些reads比对到参考基因组或者转录组。在比对中通常允许适当的错配和机构差异。目前已有的比对分析软件有（Bowtie、BWA、TopHat等）。
  • 接下来进行转录组的重建，即利用reads位置信息推断表达转录本的外显子结构，组装出转录单元。方法上主要分为两类：基因组引导法（依赖于参考基因组或转录组，如Cufflinks）、基因组独立法（也称为从头装配，是基于图论的思想，如Trinity）。
  • 接下来是对转录本的表达水平进行定量，因为在测序过程中，较长的转录本容易产生较多的reads，而且每次测序轨道上产生的reads总数也不同，因此需要对reads计数进行标准化。常用的测度是RPKM。
  • 最后进行差异表达分析，常用的软件包括DESeq、Cuffdiff、edgeR。DESeq和edgeR是基于负二项分布的统计学模型，利用Fisher精确检验方法来进行差异表达检验。Cuffdiff是基于t检验方法评估差异表达的显著性。
  1. 补充：DNA测序技术的发展历史
  • 早期测序：第一个蛋白质测序是用Sanger在20世纪50年代早期测定的胰岛素序列。20世纪60年代的RNA测序首先将RNA序列分解为RNases，接下来的片段被色谱和电泳分离，然后单个片段被连续的核酸外切酶消化破译，最后由重叠部分推导出序列
  • DNA测序的发明：1976年前后，有两种方法可以在一个下午对数百个碱基进行解码。分别是有Sanger和Coulson开发的链终止过程，以及由Maxam和Gilbert使用的化学切割过程。Sanger的方法包括用DNA聚合酶标记引物的四种延伸，每一种都有微量的链终止核苷酸，以产生不同长度的片段。
  • 扩展到人类基因组：随着“鸟枪测序法”成为了人类基因组计划的主力，人类基因组的大片段被克隆成细菌人工染色体（BACs）。个体克隆被选中并生长，DNA被分离。将纯化的的DNA作为Sanger测序的模板，从凝胶扫描图像中提取信号，并调用碱基最终生成序列。
  • 大规模并行测序：在20世纪80年代和90年代，有几个小组探索了电泳测序的替代方法。尽管在HGP完成后的10年内，“大规模并行”或“下一代”DNA测序几乎完全取代了Sanger测序。
  • 实时单分子测序：最近产生了两类实时单分子测序平台。第一种方法是Pacific Biosciences（PacBio）开发的实时观察聚合酶介导的合成。PacBio的通量仍然比最高通量的NGS平台（如Illumina）要低很多，但距离NGS平台的地位并不遥远。PicBio的最小偏差，随机误差。第二种方式是纳米孔测序。
  • 从科学历史的长远来看，DNA测序仍然是一项年轻的技术。自从发明光学显微镜以来，已有大约400年的时间，这是继续使用和发展的技术。相比之下，DNA测序发明只有40年，其技术在未来几十年和几个世纪可能还将继续发展。基于转化生物医学研究的速度，并开始转变临床医学，DNA测序的影响甚至将超过显微镜。

Perl Language

实用报表提取语言（Practical Extraction and Report Language，Perl）是一种简单易用、处理速度快的编程语言，既克服了C、C++等低级语言难以编码的缺点，又避免了高级语言速度慢的缺陷。Perl擅长处理文本，代码紧凑，代码长度一般是执行相同功能的C语言的1/4-3/4

标量数据和简单的perl语句

  • 标量数据是Perl语言最简单的数据类型
  • 主要包括数字、字符串和标量变量。
  • 标量数据可通过操作符连接得到新的标量数据
  • 标量数据的值可以存储在标量变量中。
  • 标量既可以从文件或标准输入中读取，也可以写进去。
  1. 数字操作：整数和浮点数都作为双精度浮点数处理。提供四则运算符+-*/、取余运算符%、指数操作符**
  2. 字符串操作：

What is index

rnabio.org-MiniLecture

  • Index可以有很多种不同的含义，比如，在测序过程中，它可以指unique barcodes，标记不同的DNA。而在做生信分析时，它是一种快速存取的策略，所以也会在很多不同的步骤中出现。
  • 比如，我们在做基因组分析时，测序信息的文件是相当大的，每次scan这个文件十分耗费时间，index可以构建一个内容表格，这样我们可以通过index直接跳到特定的位置。
  • 虽然计算index的时间虽然很长，但只需要做一次。当然，不同的应用可能会需要不同的index策略。
  • .fai是一种fasta文件的index文件，它分别包括了contig name(片段名称)、bases in contig(片段包含碱基数)、byte index of the file where the contig begins(片段起始位置的index字节)、bases per line(每行碱基数)、bytes per line(每行字节数)

bytes per line是什么意思？如何理解？

Example index applications and associated files

Source file Indexed file Indexing tool Use case 使用情景
.bam .bai samtools index Visualize bam in IGV 在IGV中可视化
.fasta .fai faidx Extract specific sequences form ref genome 从参考基因组中提取特定序列
.vcf vcf.gz.tbi bgzip/tabix Pull out specific variants 取出特定变体
.bed .bed.gz.tbi bgzip/tabix extract specific genomic regions 提取特定的基因组区域
  • Index在比对中也是必要的，因为找出基因组中read的位置通常是不切实际的，因此，所有的读段比对方法都使用了某种index
  • 对于一个参考基因组来说，index必须且仅需构建一次，从而每次比对的时候都可以使用。不同的比对器会使用不同的index scheme，不兼容。

对表格中的几种文件类型有进一步的了解。

How to set environment variables

使用export命令设置环境变量：

1
export RNA_HOME=~/workspace/rnaseq

如果想下次进入还使用这些变量，可以把它写进.bashrc里面：

1
nano ~/.bashrc

然后使用ctrl+o保存，Enter确认文件名不变，ctrl + x退出

1
2
# 更新环境变量
source ~/.bashrc