要求

众多具备此功能的软件同样存在，但在此恳请大家先行自行查找数篇教程，初学者应统一采用htseq-count工具，该工具将为每个样本生成一个表达量数据文件。

为了构建表达矩阵，必须编写脚本将所有样本文件进行整合。可参考生信编程课程中的第四题，即如何将多个相同的行列式文件进行合并操作。

您可以在Excel或R软件中自行尝试，进行求取平均值和计算方差等操作。

观察那些在生物学上具有特殊意义的基因表达情况，例如GAPDH和β-ACTIN等。

理论基础

在上一篇文章的对比分析中，我们面临着一个抉择：是否真的需要进行这种比对。若你仅是想了解现有基因的表达状况，那么选用-free型工具（例如，某工具）便足够了。而若你的目标是搜寻novel基因，那么你将需要借助-based型工具（比如，STAR）。然而，当我们进入本篇关于基因（转录本）的定量分析时，需要考虑的因素变得更为复杂，以至于我难以用简洁的语言清晰地阐述其中的逻辑。

定量分为三个水平

在基因分析领域，我们经常使用HTSeq-count等软件。以HTSeq-count为例，这类软件的主要任务是依据read序列及其在基因中的位置，确定该read序列的归属。需要指出的是，不同的工具在处理read序列时方法各异，比如HTSeq-count就将其视为无关紧要的信息。而则是一人一份，平均分配。

在完成基因计数并得到相应的count数值后，后续分析时必须关注标准化这一环节。若需对比同一样本中不同基因的表达水平，必须考虑转录本的长度因素；转录本越长，检测到的片段就越多，这种情况下直接比较就如同让小孩与大人进行赛跑一般，缺乏可比性。在对比不同样本中同一基因的表达水平时，即便无需关注转录本的长度，也需留意测序的深度，因为测序深度越大，检测到的可能性自然也越高。此外，还需关注由GC含量引起的偏差以及测序设备可能存在的系统误差。目前，针对read count的标准化算法包括RPKM（SE）、FPKM（PE）、TPM和TMM等，这些算法间的差异及换算方式已有文献进行汇总和讨论。然而，部分下游分析软件在输入时要求使用未经标准化的原始count数据。因此，在使用相关软件时，需留意前一步骤所使用的软件是否已对数据进行标准化处理。

在转录本分析中，我们通常依赖某款经典工具及其后续版本。这些工具面临的一大挑战是，转录本亚型之间往往存在交叉重叠。特别是在二代测序的读长短于转录本长度的情况下，如何准确区分它们？幸运的是，我们拥有了三代测序技术。这些工具普遍采用的方法是（EM）。

这些软件均以某种方式为基础，目前市面上存在众多免费的软件，例如某些，它们可以跳过比对环节，直接计算出read count，从而在运行效率上有所提升。然而，近期一篇研究文献指出，这类方法在估计丰度时可能会出现样本特异性和读长偏差的问题。

在外显子使用情况上，其统计结果与基因水平相近。然而dnastar key，值得注意的是，为了确保计数准确，我们必须提供无重叠的外显子区域的gtf文件。为此，我们提供了一款名为（ation.py）的脚本，用于执行这一分析差异外显子使用的任务。

小结

基因层面的计数、转录本层面的计数以及外显子使用层面的计数，构成了三个不同的计数等级。

标准化方法： FPKM RPKM TMM TPM

输出表达矩阵

在RNA测序分析过程中，每个基因都代表一个特征，其本质是所有外显子的总和。在可变剪接分析中，我们可以将每个外显子独立考虑。至于ChIP测序分析，则是基于预先设定的结合域。然而dnastar key，确定一个读段究竟归属何方往往颇具挑战。鉴于此，HTSeq提供了三种不同的处理模式，具体示意图请参考前一幅图。

基本用法非常的简单：

# 安装
conda install htseq
# 使用
# htseq-count [options]  
使用htseq-count工具，以正方向读取模式进行计数，针对BAM格式的文件RNA-Seq/aligned/SRR3589957_sorted.bam，对参考基因组文件reference/gencode.v26lift37.annotation.sorted.gtf进行比对，并将结果输出至SRR3589957.count文件中。

通过循环结构对位于/mnt/f/Data目录下的多个BAM文件进行操作处理。

mkdir -p RNA-Seq/matrix/
for i in `seq 56 58`
do
执行htseq-count命令，设置参数为-s no，-r pos，-f bam，读取RNA-Seq/aligned目录下SRR35899${i}_sorted.bam文件，与reference/gencode.v26lift37.annotation.sorted.gtf文件进行比对，并将结果输出到RNA-Seq/matrix目录下，生成SRR35899${i}.count文件；同时，将错误信息重定向到RNA-Seq/matrix目录下的SRR35899${i}.log文件。
done

操作过程可能较为耗时，因此不妨研究一下各个参数的具体应用，特别是以下这些使用频率较高的：

"ENSG.5_2"，"ENST.2_1"，空白字符，空白字符，空白字符，"002"，数字2，"ENSE.1_1"，二级水平，……

Jimmy的伏笔

本次研究基于人类mRNA-Seq测序数据，然而实际下载的sra文件仅有三个。通常情况下，RNA-Seq的数据分析至少需要两组重复数据，因此理论上至少需要四个sra文件。在认真阅读完关于实验方法的那一部分后，我发现其中有一部分数据是引用了其他研究小组的成果：针对293细胞，我们使用了mRNA测序技术，其数据包括（1）我们自行对-Flp-In T-REx 293细胞进行的测序结果，以及（2）由其他研究小组对-Flp-In T-REx 293细胞进行测序后提供给我们（即共享）的数据。在与Jimmy的沟通中，他坦白仅对小鼠的数据进行了研究，并未涉及人类数据。因此，我们必须依照文章中的提示，下载另一份数据资料，以便继续开展后续分析工作。

此时，人们常会提出一个疑问：能否直接对比来自不同渠道的RNA-Seq数据？更具体地，若这些数据的文库构建方法各不相同，我们又该如何进行对比呢？在进行不同来源RNA-Seq结果间的比较时，必须考虑到批次效应所带来的影响。

在应对批次效应问题时，根据我所查阅的资料，发现FPKM/RPKM这一标准化方法并不适用。对于这一标准化的批评，我在搜索过程中找到了以下看法：

有人推荐了一种名为“http://www..org//devel/bioc/html/sva.html”的软件包，我尚未深入探究其详情，希望在有限的生命里能够去学习和完善相关知识。

有人引用了文献IVT-seq中的偏差问题，来阐述RNA-Seq在不同文库间所呈现出的显著结果差异。

结论：不妨咨询一下原作者，了解他们是如何对数据进行处理的。令人惊讶的是，即便是没有重复的分析，竟然也能顺利通过审核。我们可以采用小鼠的数据来进行分析，或者选择一种不重复的分析方法，亦或是构建一份模拟数据，确保整个流程得以顺利完成。

合并表达矩阵

HTSeq-count生成的输出是单独的文件，后续处理步骤要求将这些文件整合为一个表格，其中包含行为基因名作为一列，样本名作为另一列，以及位于中间的count行列式文件。显然，手动使用Excel进行操作是不现实的（尽管可以编写一个VBA脚本，但面对庞大的数据量，这种方法并不适用），因此，此处选择编写一个脚本来完成这项任务。

基本逻辑：

读取文件

若字典中缺少对应的键，则创建新的键值对；若键已存在，则在原有值的基础上进行添加。

输出

#!/usr/bin/python3
import sys
mydict = {}

for file in sys.argv[1:]:
在遍历文件内容时，针对每一行，通过打开文件并使用读取模式，执行以下操作：
将行内容去除两端空白字符后，通过制表符分隔，得到键和值。
        if key in mydict:
在mydict字典中，将键key对应的值增加一个制表符，即执行mydict[key] = mydict[key] + '\t' + value操作。
        else:
字典中对应键的值被设置为value。
输出关键信息，同时以制表符分隔键与值，随后换行。

这段代码中设置了两个番茄钟的计时，然而调试过程耗费了一个番茄钟的时间。问题主要源于对字符串和列表这两种数据类型的混淆使用。此外，还存在一个错误：由于词典是无序的，原本标记样本来源的第一行信息，竟然错位出现在了其他行。

在论坛上搜寻到了一段更为精炼的代码，该代码遵循基因名顺序排列的要求，并成功应用了paste、awk以及shell命令。

将所有.txt文件内容输入，经过awk命令处理，首先输出每行第一个字段，然后使用制表符分隔，接着对剩余字段进行遍历，从第二个字段开始依次输出。<=NF;i+=2) printf $i"\t";printf $i}'

将文件保存为.py格式，输入源文件名为count，执行后结果将直接输出至标准输出流，并需进行输出重定向操作。

简单分析

这一步需要用到R语言或者是Excel读取数据。

1.导入数据

options(stringsAsFactors = FALSE)若样本量较小，请导入数据。control <- read.table(F盘下的Data文件夹中，RNA-Seq矩阵文件，名为SRR3589956.count的文件。,
                       sep="\t", col.names = c("gene_id","control"))
rep1 <- read.table(F盘中的Data文件夹里，RNA-Seq子目录下，matrix子目录内，存在一个名为SRR3589957的count文件。,
                    sep="\t", col.names = c("gene_id","rep1"))
rep2 <- read.table(在F盘的Data文件夹中，有一个名为RNA-Seq的子文件夹，该文件夹内又包含一个matrix子文件夹，在其中可以找到名为SRR3589958的计数矩阵文件。,
                    sep="\t",col.names = c("gene_id","rep2"))

在注释文件中，例如ENSG.13_3这样的标识在EBI数据库中无法被检索到，因此我决定仅保留ENSG这一部分信息。

合并数据，并移除无用的行。
raw_count <- merge(merge(control, rep1, by="gene_id"), rep2, by="gene_id")
raw_count_filt <- raw_count[-1:-5,]
ENSEMBL <- gsub("(.*?)\\.\\d*?_\\d", "\\1", raw_count_filt$gene_id)
row.names(raw_count_filt) <- ENSEMBL

3.总体情况, 大部分基因都为0，所以可以删掉节省体积。

summary(raw_count_filt)
严格控制，确保rep1和rep2得以有效执行。
 Min.   :     0.0   Min.   :     0.0   Min.   :     0.0  
 1st Qu.:     0.0   1st Qu.:     0.0   1st Qu.:     0.0  
 Median :     0.0   Median :     0.0   Median :     0.0  
 Mean   :   356.1   Mean   :   370.3   Mean   :   316.6  
 3rd Qu.:    15.0   3rd Qu.:    15.0   3rd Qu.:    10.0  
 Max.   :161867.0   Max.   :121902.0   Max.   :105565.0

4.看几个具体基因

在EBI上搜索GAPDH找到ID为ENSG。

GAPDH <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000111640",]
                             gene_id control  rep1  rep2
ENSG00000111640 ENSG00000111640.14_2   41857 53902 55302

文章研究的（ENSG）的表达量在KD中都降低了

> AKAP95 <- raw_count_filt[rownames(raw_count_filt)=="ENSG00000105127",]
> AKAP95
基因标识符、控制样本、重复实验1、重复实验2
ENSG00000105127，编号为ENSG00000105127.8_2，其基因序列长度为1168个碱基，其中编码区占据539个碱基，非编码区则由506个碱基组成。

下面的差异基因表达，让我想想，该如何收拾Jimmy挖的坑。

如有侵权请联系删除！