精华内容
下载资源
问答
  • 欢迎关注”生信修炼手册”!得到基因/转录本的表达量之后,通常会通过以下三种类型的图表来检验和分析生物学样本和实验设计间关系。1. 样本的聚类树利用所有样本的表达量数据,对样本进行聚类。...

    欢迎关注”生信修炼手册”!

    得到基因/转录本的表达量之后,通常会通过以下三种类型的图表来检验和分析生物学样本和实验设计间关系。

    1.  样本的聚类树

    利用所有样本的表达量数据,对样本进行聚类。理论上如果样本和实验操作都没有问题,那么属于同一组的生物学重复样本会聚到一起。示意图如下

    上图中,样本的名称用组别代替,可以看到,同一条件的样本聚在了一起。

    2. PCA图

    通过主成分分析进行降维,在二维或者三维平面上展示样本点的分布,根据点的位置,也可以看出属于同一组的样本是否在一起,不同组之间的样本有没有明显分开,示意如下


    从图中可以看到,不同条件的样本区分的很明显,而生物学重复之间距离较近,表明生物学重复的一致性和不同分组的差异性较好。

    3.  热图

    相比样本的聚类树,热图包含了更多的信息,比如可以直观的展示不同分组间表达量的差异,也是常见的可视化手段之一,示意如下

    只要有样本的表达量矩阵,DESeq2可以轻松的画出以上3种图表。但是我们应该选择原始的表达量矩阵,还是归一化之后的表达量矩阵来画呢?或者有没有其他的选择呢?

    输入的矩阵不同,得出的结论也会不同。由于基因的表达水平在不同样本间本身就存在一定的差异,所以无论是采用原始的还是归一化之后的表达量矩阵,效果都不理想。针对这一问题,DESeq2提出了两种count值的转换算法,rlogVST转换。

    1. rlog 转换

    rlog 转换的用法如下

    rld <- rlog(dds)

    2. VST 转换

    用法如下

    vsd <- vst(dds)

    两种转换本质上是在降低生物学重复之间的差异,使得样本聚类和PCA分析的效果更好。转换之后的表达量数据可以采用assay函数进行提取,代码如下

    > head(assay(rld)[, 1:2])
           sample1   sample2
    gene1 2.049029 1.6828707
    gene2 8.151262 6.8552583
    gene3 0.818971 0.2964686
    gene4 5.340361 4.4766682
    gene5 6.316175 6.8345783
    gene6 2.157821 1.9264385

    对于raw count定量表格,建议采用rlog或者VST转换之后的数据去进行PCA和聚类分析,效果会更好。

    利用DESeq2提供的示例数据pasilla,分别用原始的count, 归一化之后的count, rlog, vst 转换的count 进行PCA分析,代码如下

    dds <- estimateSizeFactors(dds)
    raw <- SummarizedExperiment(counts(dds, normalized=FALSE),
                                    colData=colData(dds))
    nor <- SummarizedExperiment(counts(dds, normalized=TRUE),
                                    colData=colData(dds))
    vsd <- vst(dds)
    rld <- rlog(dds)
    pdf("PCA.pdf")
    plotPCA( DESeqTransform(raw), intgroup=c("condition", "type") )
    plotPCA( DESeqTransform(nor), intgroup=c("condition", "type") )
    plotPCA(vsd, intgroup=c("condition", "type"))
    plotPCA(rld, intgroup=c("condition", "type"))
    dev.off()

    raw count 的结果如下


    归一化之后count结果如下

    VST转换之后的结果如下

    rlog转换之后的结果如下

    可以很明显看出,原始的count和归一化之后的count, 其PCA图是杂乱无序的,没什么明显规律,而VST和rlog转换之后,生物学重复之间更佳的接近,不同分组也区分的较为明显。

    ·end·

    —如果喜欢,快分享给你的朋友们吧—

    扫描关注微信号,更多精彩内容等着你!

    展开全文
  • 当转录组研究的样本分组超过2组时,...表达模式聚类热图以各样本中基因的表达量绘制热图,在图中每列表示一个样本,每行表示一个基因,图中的颜色的深浅表示基因在该样本中的表达量。表达模式聚类热图图像左侧的聚类...

    5db89f626f57d59f3ac4fde1c4afc7c6.png

    当转录组研究的样本分组超过2组时,由于差异表达基因的识别只能针对两组样本,因此,只通过差异表达基因无法分析整体上的基因表达变化规律,此时可以通过表达模式聚类将基因按照其在不同样本中的表达变化规律进行归类,进而推测其与特定功能的可能联系

    表达模式聚类热图

    以各样本中基因的表达量绘制热图,在图中每列表示一个样本,每行表示一个基因,图中的颜色的深浅表示基因在该样本中的表达量。

    6d02efb0b6ef11d7cc3ca3f87ed859b7.png
    表达模式聚类热图

    图像左侧的聚类数表示各个基因在所有样品中表达规律的相似性,在聚类中分支越近的基因,其表达量的变化规律就越接近,同一聚类模式的基因可能具有相同或相关的功能

    表达模式聚类折线图

    表达模式聚类折线图的绘制过程如下:首先根据基因在不同样本中的表达量将其分为多个subcluster,之后每个subcluster分别进行绘图,图像中x轴代表不同的样本,y轴为各个基因表达量的对数值。

    dcd4f8b9f2110c7ebcc9f5582df3d80f.png
    表达模式聚类折线图

    灰色线条表示一个subcluster中的基因在不同样本中的相对表达量 (聚为一类的基因会有很多个,图中显示比较密集)。

    蓝色线条表示这个subcluster中的所有基因在不同样本中相对表达量的平均值。

    表达模式聚类热图通常是按照欧式距离等相关性距离计算方法进行聚类,而折线图可以根据需要人为指定分类的数目。

    Venn分析

    分别识别不同组样本两两之间的差异表达基因,之后进行Venn分析。

    44c64a5727cedcd424012500bdc96e4d.png
    差异基因Venn分析

    在图中,不同颜色的圆圈代表两组样本的差异基因,圆内部所有数字之和代表该组比较中差异基因个数的总和,圆的交叉区域代表不同比较的差异基因中共有基因个数,而最外圈单独的部分表示该组比较中特有的差异基因数目。

    该分析可以识别在研究的所有过程中均发挥作用的关键基因 (所有比对组中的共有基因),也可用于识别与特定功能相关的关键基因 (各组比较中的特有基因)

    a2e56bc2e47d2bfe39c96c6e144f0cf4.png
    展开全文
  • R在读取和处理数据的过程中会将所有的变量和占用都储存在RAM当中,这样...但是最近我发现了一个基于python的单细胞基因表达分析包scanpy,能够很好地在我这个仅4G内存的小破机上分析378k的细胞,并且功能丰富程度...

    R在读取和处理数据的过程中会将所有的变量和占用都储存在RAM当中,这样一来,对于海量的单细胞RNA-seq数据(尤其是超过250k的细胞量),即使在服务器当中运行,Seurat、metacell、monocle这一类的R包的使用还是会产生内存不足的问题。

    但是最近我发现了一个基于python的单细胞基因表达分析包scanpy,能够很好地在我这个仅4G内存的小破机上分析378k的细胞,并且功能丰富程度不亚于Seurat。它包含了数据预处理、可视化、聚类、伪时间分析和轨迹推断、差异表达分析。根据官网描述,scanpy能够有效处理超过1,000,000个细胞的数据集。

    Scanpy – Single-Cell Analysis in Python:https://scanpy.readthedocs.io/en/latest/

    安装与数据下载

    安装好python3之后,终端运行:

    pip install scanpy

    Step0, 读取数据

    运行python

    import numpy as np

    import pandas as pd

    import scanpy as sc

    # 可以直接读取10Xgenomics的.h5格式数据

    adata = sc.read_10x_h5("/Users/shinianyike/Desktop/ica_bone_marrow_h5.h5"

    , genome=None, gex_only=True)

    adata.var_names_make_unique()

    查看数据:

    adata

    AnnData object with n_obs × n_vars = 378000 × 33694

    var: 'gene_ids'

    共378000个细胞,33694个基因。

    Step1, 数据预处理

    这一步目的将数据进行筛选和过滤,一些测序质量差的数据会让后续的分析产生噪音和干扰,因此我们需要将它们去除。

    展示在所有的细胞当中表达占比最高的20个基因。

    sc.pl.highest_expr_genes(adata, n_top=20)

    表达水平前20的基因.png

    基础过滤:

    去除表达基因200以下的细胞;

    去除在3个细胞以下表达的基因。

    sc.pp.filter_cells(adata, min_genes=200) # 去除表达基因200以下的细胞

    sc.pp.filter_genes(adata, min_cells=3) # 去除在3个细胞以下表达的基因

    adata

    AnnData object with n_obs × n_vars = 335618 × 24888

    obs: 'n_genes'

    var: 'gene_ids', 'n_cells'

    共335618个细胞,24888个基因。(可以发现细胞跟基因数量都减少了)

    质量控制:

    在测序后,有很大比例是低质量的细胞,可能是因为细胞破碎造成的胞质RNA流失,由于线粒体比单个的转录分子要大得多,不容易在破碎的细胞膜中泄漏出去,这样一来就导致测序结果显示线粒体基因的比例在细胞内占比异常高。这一步质量控制的目的就是将这些低质量的细胞去除掉。

    计算线粒体基因占所有基因的比例:

    mito_genes = adata.var_names.str.startswith('MT-')

    # for each cell compute fraction of counts in mito genes vs. all genes

    # the `.A1` is only necessary as X is sparse (to transform to a dense array after summing)

    adata.obs['percent_mito'] = np.sum(

    adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1

    # add the total counts per cell as observations-annotation to adata

    adata.obs['n_counts'] = adata.X.sum(axis=1).A1

    作小提琴图,查看线粒体基因占比分布:

    sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'],

    jitter=0.4, multi_panel=True)

    细胞表达的基因种数、基因数、线粒体基因占比

    由于数据点实在太多,小提琴已被数据点覆盖,无法显示出来。

    这里还可以做一个散点图,也可以直观地显示出一些异常分布的数据点。

    sc.pl.scatter(adata, x='n_counts', y='percent_mito')

    sc.pl.scatter(adata, x='n_counts', y='n_genes')

    adata

    AnnData object with n_obs × n_vars = 335618 × 24888

    obs: 'n_genes', 'percent_mito', 'n_counts'

    var: 'gene_ids', 'n_cells'

    335618个细胞,24888个基因;

    下面这一步进行真正的过滤,根据上面做的图,去除基因数目过多(大于等于4000)和线粒体基因占比过多(大于等于0.3)的细胞:

    adata = adata[adata.obs['n_genes'] < 4000, :]

    adata = adata[adata.obs['percent_mito'] < 0.3, :]

    过滤后查看剩下多少细胞和基因。

    adata

    View of AnnData object with n_obs × n_vars = 328435 × 24888

    obs: 'n_genes', 'percent_mito', 'n_counts'

    var: 'gene_ids', 'n_cells'

    328435个细胞,24888个基因。

    数据标准化

    在绘图之前,还要进行一步数据标准化,将表达量用对数计算一遍,这样有利于绘图和展示。

    sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)

    sc.pp.log1p(adata)

    adata.raw = adata # 储存标准化后的AnnaData Object

    识别差异表达基因

    sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

    sc.pl.highly_variable_genes(adata)

    将保守的基因去除,留下差异表达的基因用于后续分析。

    adata = adata[:, adata.var['highly_variable']]

    adata

    View of AnnData object with n_obs × n_vars = 328435 × 1372

    obs: 'n_genes', 'percent_mito', 'n_counts'

    var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'

    328435个细胞,1372个基因。

    回归每个细胞总计数和线粒体基因表达百分比的影响。将数据放缩到方差为1。单细胞数据集可能包含“不感兴趣”的变异来源。这不仅包括技术噪音,还包括批次效应,甚至包括生物变异来源(细胞周期阶段)。正如(Buettner, et al NBT,2015)中所建议的那样,从分析中回归这些信号可以改善下游维数减少和聚类。

    这一步高内存需求预警,建议清理电脑缓存,关闭后台不使用的应用。

    sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])

    /Users/shinianyike/anaconda3/lib/python3.6/site-packages/statsmodels/compat/pandas.py:56: FutureWarning: The pandas.core.datetools module is deprecated and will be removed in a future version. Please use the pandas.tseries module instead.

    from pandas.core import datetools

    Scale each gene to unit variance. Clip values exceeding standard deviation 10.

    sc.pp.scale(adata, max_value=10)

    Step2, 主成分分析

    主成分分析是一种将数据降维的分析方法,是考察多个变量间相关性一种多元统计方法,研究如何通过少数几个主成分来揭示多个变量间的内部结构,即从原始变量中导出少数几个主成分,使它们尽可能多地保留原始变量的信息,且彼此间互不相关.通常数学上的处理就是将原来P个指标作线性组合,作为新的综合指标。

    sc.tl.pca(adata, svd_solver='arpack') # PCA分析

    sc.pl.pca(adata, color='CST3') #绘图

    作碎石图,观测主成分的质量。这个图用于选择后续应该使用多少个PC,用于计算细胞间的相邻距离。

    sc.pl.pca_variance_ratio(adata, log=True)

    在这里先将数据保存,便于后续操作的更改。

    adata.write("pca_results.h5ad")

    Step3, 聚类分析

    计算细胞间的距离:

    这里的参数就先按照默认值设定:

    sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

    参数说明:

    n_neighbors指的是每个点的邻近点的数量,据评论区@小光amateur 所说neighbors的个数越多,聚类数会越少。

    pc的数量依赖于上面所做的碎石图,一般是选在拐点处的的主成分,只需要一个粗略值,不同的pc数量所产生的聚类形状也不同。我后来更改为PC=16,效果比下图要好一些。

    将距离嵌入图中:

    sc.tl.umap(adata)

    sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

    sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)

    聚类

    sc.tl.louvain(adata)

    sc.pl.umap(adata, color=['louvain'])

    output_51_0.png

    这里得到了29类细胞,每个颜色代表一种。

    将数据保存。

    adata.write("umap.h5ad")

    寻找marker基因

    marker基因通常是细胞表面抗原,用于定义出该细胞的细胞类型。

    为了定义每个簇属于什么细胞,根据基因的差异表达水平,将每个簇排名前25的基因导出。

    sc.tl.rank_genes_groups(adata, 'louvain', method='wilcoxon')

    sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

    下一步的工作是找出每一个簇的marker基因对应的细胞类型,这主要依靠一些数据库或生物学的相关背景知识。

    文章已发布到微信公众号:百味科研芝士,欢迎关注。

    展开全文
  • 基因表达模式聚类以及可视化

    万次阅读 2017-07-19 21:44:18
    最近在使用RNA_seq数据做些分析,结果得到了大量差异表达以及共表达的基因,如何合理展示这些基因...输入文件就是一个基因表达量的矩阵,如下图。 代码也很简单,见下图,也请点击阅读原文查看代码。​ library(TCseq)

    最近在使用RNA_seq数据做些分析,结果得到了大量差异表达以及共表达的基因,如何合理展示这些基因也是一件不简单的事情。除了常见的热图(heatmap)展现形式,今天在推荐另外一种展示方式(下图C)。

    这里写图片描述

    需要R包TCseq或者Mfuzz。我这里给出的代码是基于TCseq。
    输入文件就是一个基因表达量的矩阵,如下图。
    这里写图片描述

    代码也很简单,见下图,也请点击阅读原文查看代码。

    library(TCseq)
    svg(filename="/stress_analysis/drought_and_heat/heat.svg",width=3.5,height=15)
    data <- read.csv("/stress_analysis/drought_and_heat/heat.txt",row.names="Name",sep="\t")
    data <- as.matrix(data)
    tca <- timeclust(data, algo = "cm", k = 6, standardize = TRUE)
    p <- timeclustplot(tca, categories="time(h)", cols = 1, axis.text.size=12,legend.text.size=10, legend.title.size=10,title.size=12,axis.title.size=12)
    cluster1 <- clustCluster(tca) #输出基因的聚类信息
    out_file <- paste("//stress_analysis/drought_and_heat/heat_cluster.txt",sep ="")
    write.csv(cluster1,out_file)
    dev.off()

    结果如下图,该图经过后期编辑,与热图进行了合并。

    这里写图片描述

    展开全文
  • 聚类、分类、回归、关联分析 1 classification (分类)。 分类是找出数据库中的一组数据对象的共同特点并按照分类模式将其划分为不同的类,其目的是通过分类模型,将数据库中的数据项映射到摸个给定的类别中。可以...
  • 回归分析反映了数据库中数据的属性值的特性,通过函数表达数据映射的关系来发现属性值之间的依赖关系。它可以应用到对数据序列的预测及相 关关系的研究中去。在市场营销中,回归分析可以被应用到各个方面
  • 例如,数以百万计的车站、街道、机场和世界各地的城市的摄像头所产生的图像数据、海量的股票信息、大规模的高维度基因表达数据、商品信息、文档信息等,这对于压缩、存储、聚类和传输大量复杂的高维数据提供了原材料...
  • 回归分析反映了数据库中数据的属性值的特性,通过函数表达数据映射的关系来发现属性值之间的依赖关系。它可以应用到对数据序列的预测及相关关系的研究中去。在市场营销中,回归分析可以被应用到各个方面。如通过
  • K-均值聚类

    千次阅读 2010-06-11 11:52:00
    该算法原理简单并便于处理大量数据,在基因表达数据分析中 得到广泛应用,如Tavazoie等应用K-means聚类酵母细胞周期表达数据。在K-means算法运行前必须先指定聚类数目K和迭代次数或收敛条 件,并指定K个初始聚类...
  • 热图(Heatmap):用颜色变化直观的表达数据之间差异的图,是对实验数据进行质制和差异数据的展现,是数据挖掘类文章的标配。Rplot1.jpeg例如上图,...上方树形图表示对来自不同实验分组的不同样品的聚类分析结果Draw...
  • 该算法原理简单并便于处理大量数据,在基因表达数据分析中 得到广泛应用,如Tavazoie等应用K-means聚类酵母细胞周期表达数据。在K-means算法运行前必须先指定聚类数目K和迭代次数或收敛条 件,并指定K个初始聚类...
  • 热图(heatmap)的典型应用是简单地聚合大量数据,并使用一种渐进的色带来优雅地表现,最终效果一般优于离散点的直接显示,可以很直观地展现空间...如果对表达量去一下log10,发现10000变成了4,10变成了1,这样之前离.
  • 该算法原理简单并便于处理大量数据,在基因表达数据分析中得到广泛应用,如Tavazoie等应用K-means聚类酵母细胞周期表达数据。在K-means算法运行前必须先指定聚类数目K和迭代次数或收敛条件,并指定K个初始聚类中心,...
  • 以RNA芯片表达数据为例来讲,即使经过一系列分析后缩小了范围后,仍旧存在多个基因表达数据需要展示,这就使得热图在展示多个表达矩阵数据时的具有很强大的数据可视化作用。热图具有大致以下几方面的优点:1)可以...
  • 圭 非负矩阵分解及其在基因表达数据分析中的应用 曹胜玉 刘来福 索藩蕊大学数学秘警学茨察 摘要介绍非负矩阵分解的基本原理及其在生物信息学中基因表达数据分析巾的应用并将该方法用子一缀白血瘸微阵翔数据的聚类黎...
  • WGCNA如何挖掘潜在的共表达基因

    千次阅读 2018-11-05 11:34:37
    共表达基因指的是表达量具有协同变化趋势的基因集合,通常认为这些基因参与相同的生物学过程,比如参与同一个代谢通路,正是由于功能上的协同作用,导致表达量呈现出高度相关性。 在WGCNA中,对传统的相关系数进行...
  • 基因表达数据中信息基因和基因调控网络 第六周报告 本周主要看了《基因芯片技术》《基因表达数据的聚类分析》两篇论文,初步了解了基因芯片和聚类分析的含义。 一、基因芯片技术基因芯片技术是同时将大量的探针分子...
  • 趋势分析--转载

    2017-09-11 15:04:00
    这个软件的主要用途就是针对时间顺序取样的表达量数据进行聚类分析其表达模式。如果你使用我们的趋势分析工具完成分析,在撰写文章的时候可以直接引用这款STEM软件[1]。 我们都了解RNA-seq的基础分析...
  • 趋势分析:STEM软件

    2021-05-10 15:45:27
    这个软件的主要用途就是针对时间顺序取样的表达量数据进行聚类分析其表达模式。如果你使用我们的趋势分析工具完成分析,在撰写文章的时候可以直接引用这款STEM软件[1]。 我们都了解RNA-seq的基础分析是差异表达。...
  • Origin 如何做主成分分析

    千次阅读 2020-11-01 19:57:01
    在高通量测序中,主要基于基因表达量、种群丰度等进行样本的聚类,下图是一篇客户文章的基于表达量的PCA结果。 (Plant biotechnology journal, 2018) 那么该如何进行主成分分析呢?今天为大家介绍如何用...
  • 基因芯片技术,也称为核酸阵列芯片技术,产生于20世纪90年代,在近...通过有效的基因搜索和相关基因表达谱的聚类,提取表达显著差异的基因和共表达基因; pd=gprread(‘mouse_h3pd.gpr’) figure maimage(pd,‘F635 Med
  • 应用主成分分析和层次聚类方法。 所需的包 安装 您可以使用 install_github()函数从安装当前版本的PPClustA: if ( ! requireNamespace( " devtools " , quietly = TRUE )) install.packages( " devtools " ) ...
  • Q&A (第三期)

    2020-12-09 11:11:15
    何时采用聚类分析? 考察特征 样本很大的时候考虑什么聚类? k均值 最近距离法、最远距离法、重心法和类平均法的英文表达 略 Kmeans聚类的R实现函数为? kmeans 如何说明每一类的特征或特点? kmeans给出重心 连线...
  • matlab神经网络30个案例分析

    千次下载 热门讨论 2011-06-01 20:06:07
    根据货运影响因素的分析,分别取国内生产总值(GDP),工业总产值,铁路运输线路长度,复线里程比重,公路运输线路长度,等级公路比重,铁路货车数量和民用载货汽车数量8项指标因素作为网络输入,以货运总量,铁路...
  • 基于多尺度几何分析方法——非下采样轮廓波(Contourlet)变换(NSCT)和Beamlet变换,提出...数值实验表明,与传统的融合方法相比较,本文方法能够有效减少噪声对融合图像的干扰,增强了融合的线性细节表达能力,提高了信息
  • 为了在现有的手势数据基础上合成新的地方手语手势,采取对这些手势数据进行运动相似性分段和动态聚类自动获取具有时序特征的基本手形类数据以及结合手臂状态数据的方法,最后生成新的地方手语手势动画数据。...
  • 如社交网络、知识图谱、生物基因网络等,随着网络数据规模的不断增加,如何简化表达大规模网络图结构已成为图可视化领域中的研究热点,经典的网络图简化可视化方法主要包括图采样、边绑定和图聚类等技术,在减少大量...
  • ApobecEnrich_define.R:运行tSNE和期望最大化聚类的迭代,以创建Apobec突变签名富集标签,并将其与已建立的Apobec富集得分进行比较。 RandomForest_100sims.R:从平衡的APOBEC富集/非富集样本中创建100个随机森林...
  • 在上文我们已经基本完成了聚类分析和细胞在不同条件下的比例分析,下面继续开始啦,让你看看iCellR的强大! 1. 每个cluster的平均表达量 my.obj <- clust.avg.exp(my.obj) head(my.obj@clust.avg) # gene ...

空空如也

空空如也

1 2 3
收藏数 53
精华内容 21
关键字:

表达量聚类分析