合并成表达矩阵

2023-05-16

3.合并成表达矩阵
rsem官方提供了将多个样本合并成一个表达矩阵的工具，使用如下：

rsem-generate-data-matrix ./C01/C01.genes.results ./C02/C02.genes.results ./C03/C03.genes.results >output_matrix.txt

笔者并不喜欢官方提供的脚本，并且只有count矩阵结果，没有TPM和FPKM矩阵，所以可以写一个简单的python脚本搭配shell脚本重新将所需要的结果提取出来。

vim rsem-count-extract.py
#!/usr/bin/env python
#coding: utf-8
import sys
from itertools import islice
mydict = {}
for a in sys.argv[1:]:
c = open(a, ‘r’)
for line in islice(c, 1, None):
a = line.strip().split()
key = a[0]
value = a[4]
#可以选择其他列
if key in mydict:
mydict[key] = mydict[key] + ‘\t’ + value
else:
mydict[key] = value
for key,value in mydict.items():
print(key + ‘\t’ + value)
然后将python脚本的调用整合在shell脚本里，

list_csv=""
list_name=“gene”
for i in ./*/genes.results
do
echo i l i s t c s v = i list_csv= ilistcsv={list_csv}" " i i n a m e 1 = {i} i_name1= iiname1={i%/}
i_name2=KaTeX parse error: Expected '}', got '#' at position 9: {i_name1#̲#*/} list_name={list_name}" "${i_name2}
echo ${list_name}
done
echo $list_name > ./gene-count-matrix.txt
python ./rsem-count-extract.py ${list_csv} >> ./gene-count-matrix.txt
这样就可以批量提取到样本名称正确的count矩阵，同时如果想求TPM或者FPKM的矩阵，可以修改python脚本中value = a[4]的数值为5或者6。

得到表达矩阵后，后续分析可以采用官方推荐的ebseq或者比较常用的deseq2做差异基因分析。对于转录本的差异分析目前来说并不建议，这个主要还是因为现有算法对转录本定量的精确度并不高。

。。。。。。。。。。。。。。

完成之后得到这些文件，其中，rsem.genes.results和rsem.isoforms.results分别表示gene水平和转录本水平的定量结果。每一列含义：

less rsem.genes.results|head -n 1
gene_id transcript_id(s) length effective_length expected_count TPM FPKM
less rsem.isoforms.results|head -n 1
transcript_id gene_id length effective_length expected_count TPM FPKM IsoPct

后面用EBseq检验差异基因/转录本时，会使用到这两个文件。

Differential Expression Analysis using EBSeq
下面以gene水平差异分析为例。

3.1 generate-data-matrix
这一步提取上一步得到的每个样本定量结果文件中的expected_count列，组成数据矩阵。https://www.jianshu.com/p/3e96d0f8fe71

~/software/rsem/rsem-generate-data-matrix
SRR1.rsem.genes.results SRR2.rsem.genes.results
SRR3.rsem.genes.results SRR4.rsem.genes.results
SRR5.rsem.genes.results SRR6.rsem.genes.results
SRR7.rsem.genes.results SRR8.rsem.genes.results \

~/project/RNA-seq/count/GeneMat.txt

。。。。。。。。。。。。。。。。。。。

差异表达

mkdir 06deseq_out

#构建表达矩阵
rsem-generate-data-matrix *rsem.genes.results > output.matrix 转录组分析流程——STAR+RSEM+Deseq2

BeautifulSoulpy
https://www.jianshu.com/p/ad605d4fa6f6

表达矩阵的合并
https://www.jianshu.com/p/853f3706eaa8

。。。。
htseq合并计数https://www.jianshu.com/p/688a42dc07b8

。。。。。。。。.。。。.
后期分析

2020-01-23 RSEM转录组RNA-seq测序分析实操
本文涉及到的软件
RSEM软件包 RSEM (RNA-Seq by Expectation-Maximization)

Trinity软件包

安装：
conda install -c bioconda trinity

rsem-calculate-expression --strandedness reverse
–paired-end
-p 16
–bowtie2
…/cleandata/W2BR_P_R1.fq.gz
…/cleandata/W2BR_P_R2.fq.gz
/public/home/huangshsh/zxg/rat_trans_nn/rsem/ref/rn6
/path to output/W2BR

#–strandedness reverse 建库方式，一般为illumina的方式，具体自己查资料
#–paired-end 双尾测序

-p 16程序运行的线程数

–bowtie2 用bowtie2构建的索引

#…/cleandata/W2BR_P_R1.fq.gz 第一个测序文件
#…/cleandata/W2BR_P_R2.fq.gz 第二个测序文件

/public/home/huangshsh/zxg/rat_trans_nn/rsem/ref/rn6 参考基因组文件位置…/ref/,参考基因组索引前缀rn6

#W2BR 结果文件前缀（这个参数很容易遗漏，会报错的）

一个文件夹和两个文件，其中genes.results是以基因为基础的，包含了gene_id，对应transcript_id, count,TPM和FPKM等信息。
isoforms.results是转录本为基础的

正文3.生成表达量矩阵文件
利用trinity安装包中的util/abundance_estimates_to_matrix.pl（一个perl脚本）

ls *.genes.results > genes.quant_files.txt #将所有的genes.results文件的名字存到genes.quant_files.txt ，当成一个目录后边要用

/路径到trinity安装包下/util/abundance_estimates_to_matrix.pl
–est_method RSEM
–cross_sample_norm TMM
–gene_trans_map none
–quant_files genes.quant_files.txt
–out_prefix rsem_matrix

or /public/home/huangshsh/zxg/software/trinityrnaseq-2.9.0/util/abundance_estimates_to_matrix.pl --est_method --quant_files file.listing_target_files.txt

Note, if only a single input file is given, it’s expected to contain the paths to all the target abundance estimation files.

Required:

–est_method RSEM|eXpress|kallisto|salmon (needs to know what format to expect)

–gene_trans_map the gene-to-transcript mapping file. (if you don’t want gene estimates, indicate ‘none’.

可以从gtf中提取转录本信息，最后做成的文件要转录本-基因一一对应，转录本一定要在前

边，两者用\t分隔。

Options:

–cross_sample_norm TMM|UpperQuartile|none (default: TMM)

–name_sample_by_basedir name sample column by dirname instead of filename

–basedir_index default(-2)

–out_prefix default: value for --est_method 前缀随意命名，默认为RSEM

–quant_files file containing a list of all the target files.

结果文件如下图
我们后续要用到第一个，rsem_matrix.gene.counts.matrix

正文4.差异表达基因分析
要用到R包里边的DESeq2或edgeR，所以R里边要提前装好DESeq2包或者edgeR包。
脚本命令用到trinity软件目录下的/Analysis/DifferentialExpression/run_DE_analysis.pl，所用到的脚本参数如下：

$TRINITY_HOME=~/zxg/software/trinityrnaseq-2.9.0 #trinity的安装路径

perl $TRINITY_HOME/Analysis/DifferentialExpression/run_DE_analysis.pl
–matrix ~/zxg/rat_trans_nn/rsem/rsem_matrix.gene.counts.matrix
–method DESeq2
–samples_file sample.list.txt

##sample.list.txt格式如上图所示

https://blog.csdn.net/weixin_44649331/article/details/103837169?utm_medium=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.channel_param&depth_1-utm_source=distribute.pc_relevant.none-task-blog-BlogCommendFromMachineLearnPai2-3.channel_param

https://blog.csdn.net/xiaomotong123/article/details/106900481/?biz_id=102&utm_term=DESeq2%E7%9A%84design&utm_medium=distribute.pc_search_result.none-task-blog-2_allsobaiduweb~default-0-106900481&spm=1018.2118.3001.4187

https://www.jianshu.com/p/3a0e1e3e41d0

http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

condition type

treated1fb treated single-read

treated2fb treated paired-end

treated3fb treated paired-end

untreated1fb untreated single-read

untreated2fb untreated single-read

untreated3fb untreated paired-end

untreated4fb untreated paired-end

dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition) #_{在R里面用于构建公式对象，}左边为因变量，右边为自变量。
dds <- DESeq(dds) #标准化
res <- results(dds, contrast=c(“condition”,“treated”,“control”)) #差异分析结果

构建dds矩阵需要：

表达矩阵即上述代码中的countData，就是我们前面通过read count计算后并融合生成的矩阵，行为各个基因，列为各个样品，中间为计算reads或者fragment得到的整数。我们后面要用的是这个文件（mouse_all_count.txt）

样品信息矩阵即上述代码中的colData，它的类型是一个dataframe（数据框），第一列是样品名称，第二列是样品的处理情况（对照还是处理等），即condition，condition的类型是一个factor。这些信息可以从mouse_all_count.txt中导出，也可以自己单独建立。

差异比较矩阵即上述代码中的design。差异比较矩阵就是告诉差异分析函数是要从要分析哪些变量间的差异，简单说就是说明哪些是对照哪些是处理。

正式构建dds矩阵

library(DESeq2) #加载包
countData <- mouse_all_count
condition <- factor(c(“control”,“control”,“KD”,“KD”))
colData <- data.frame(row.names=colnames(countData), condition)
我们可以看一下condition和colData的具体内容：https://zhuanlan.zhihu.com/p/30350531
#正式构建dds矩阵
dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), design= ~ condition )
head(dds) #查看一下构建好的矩阵
可以看到rownames都是原来mouse_all_count.txt中的行号，我们也可以在构建dds矩阵的时候，把行号直接替换成gene_id，这样操作后生成的dds矩阵如下：

RNA-seq练习第三部分（DEseq2筛选差异表达基因,可视化）https://www.jianshu.com/p/82787ddada38

Transcript abundance files and tximport / tximeta
Our recommended pipeline for DESeq2 is to use fast transcript abundance quantifiers upstream of DESeq2, and then to create gene-level count matrices for use with DESeq2 by importing the quantification data using tximport (Soneson, Love, and Robinson 2015). This workflow allows users to import transcript abundance estimates from a variety of external software, including the following methods:

Salmon (Patro et al. 2017)
Sailfish (Patro, Mount, and Kingsford 2014)
kallisto (Bray et al. 2016)
RSEM (Li and Dewey 2011)

。。。。。。。。。。。。。。。。。。。

下面是HTseq流程的定量

转入组入门七（mac 版）：差异基因分析【把所有样本都合并成一个表如样本名称：A.rep1 A.rep2 B.rep1 B.rep2】https://www.jianshu.com/p/dc178908a4c8

读取表达矩阵

首先将四个文件分别赋值：control1，control2，rep1，rep2

control1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589959.count", sep="\t", col.names = c(“gene_id”,“control1”))
control2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589961.count", sep="\t", col.names = c(“gene_id”,“control2”))
rep1 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589960.count", sep="\t", col.names = c(“gene_id”,“akap951”))
rep2 <- read.table("/Users/chengkai/Desktop/zhuanlu_files/RNA-Seq/matrix/SRR3589962.count", sep="\t",col.names = c(“gene_id”,“akap952”))

将四个矩阵按照gene_id进行合并，并赋值给raw_count

raw_count <- merge(merge(control1, control2, by=“gene_id”), merge(rep1,rep2, by=“gene_id”))

需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter

raw_count_filt <- raw_count[-48823:-48825, ]
raw_count_filter <- raw_count_filt[-1:-2, ]

因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。

第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL

ENSEMBL <- gsub("\.\d*", “”, raw_count_filter$gene_id)
row.names(raw_count_filter) <- ENSEMBL
raw_count_filter <- raw_count_filter[ ,-1]

. 合并表达矩阵https://www.jianshu.com/p/ff585b72f04e
使用R将这四个文件进行合并（59 和 61 是对照， 60和62是处理），得到最后的融合表达矩阵，R语言代码如下

options(stringsAsFactors = FALSE)

首先将四个文件分别赋值：control1，control2，rep1，rep2

将四个矩阵按照gene_id进行合并，并赋值给raw_count

raw_count <- merge(merge(control1, control2, by=“gene_id”), merge(rep1,rep2, by=“gene_id”))

需要将合并的raw_count进行过滤处理，里面有5行需要删除的行，在我们的小鼠的表达矩阵里面，是1,2,48823,48824,48825这5行。并重新赋值给raw_count_filter

raw_count_filt <- raw_count[-48823:-48825, ]
raw_count_filter <- raw_count_filt[-1:-2, ]

因为我们无法在EBI数据库上直接搜索找到ENSMUSG00000024045.5这样的基因，只能是ENSMUSG00000024045的整数，没有小数点，所以需要进一步替换为整数的形式。

第一步将匹配到的.以及后面的数字连续匹配并替换为空，并赋值给ENSEMBL

ENSEMBL <- gsub("\.\d*", “”, raw_count_filter$gene_id)

将ENSEMBL重新添加到raw_count_filter矩阵

row.names(raw_count_filter) <- ENSEMBL

看一些基因的表达情况，在UniProt数据库找到AKAP95的id，并从矩阵中找到访问，并赋值给AKAP95变量

AKAP95 <- raw_count_filter[rownames(raw_count_filter)==“ENSMUSG00000024045”,]

查看AKAP95

AKAP95
7. 简单分析
合并后的数据框

一些名词
sample：样本（包含了实验条件和重复），例如，我们说测10个样，也就是1个处理一个对照，各5个重复
normalization:标准化，将各个样本的结果放在同一维度上，赋予它们可比性
library size normalization:文库大小标准化，用来矫正不同实验条件下的测序深度（覆盖度）。例如，条件A相对B加了2倍的材料进行测序，这时就需要平衡两者
effective length：有效长度，其实也就是长度越长的转录本，reads落在上面的数量越多，因此计算时需要平衡
gene level：在基因层次上，将每个基因视作独立的转录本，而且均包含基因的所有外显子。但是还有一些情况需要更细致的分析，因此基因表达量需要建立在转录本层次上（transcript level）
TMM(edegR) normalization: Trimmed mean of M values。edgeR软件使用的标准化算法叫TMM，它排除了一对实验条件下reads counts比率比较极端或者表达量的平均值比较极端的基因，从而估算的测序深度
DESeq normalization：DEseq包使用的算法。它选出所有基因reads counts中位数的基因，估算它的测序深度

重回正题~Trinity确实是个非常好的软件，官方文档写的十分详细，而且主流的无参转录组分析都推荐使用这款，它可以预测出更长的亚型、更多的基因和转录本，另外它还有几十个小工具，比如不会做热图，不会做火山图都可以借助它来进行可视化，完全不需要本人需要会R知识，直接一个shell命令搞定。去年2月我刚开始接触生信不久，就做了这个内容，因此记得特别清楚那种感觉——不明不白整出来一个自己看不懂的pdf，导出以后，哇这么高大上！问题：转录本ID转换？
https://www.jianshu.com/p/28b435efe79b

如果只是关注某一条基因在样本间的差异，那么看P value就好了；如果关注整体的差异情况，一般就要看FDR、q-value https://www.jianshu.com/p/b2ce13121c6c

/转录组样本间表达标准化https://www.jianshu.com/p/4ae683b257a0

转录组差异基因筛选https://www.jianshu.com/p/cebd47e84b1d

转录组差异分析金标准-Limma-voom实战https://www.jianshu.com/p/dee7346482e5

本文内容由网友自发贡献，版权归原作者所有，本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容，请联系:hwhale#tublm.com(使用前将#替换为@)