玉米中的元基因调控网络突出了功能上相关的调控相互作用。\ca.19a5.R

这段代码处理RNA-Seq数据，主要包括质量控制（QC）结果的读取、数据过滤、样本筛选、数据转换、聚类分析和降维（t-SNE）可视化。我们可以将代码分为几个部分进行讲解。

第一部分：创建工作目录并读取相关数据

yid = 'ca19a5'
dirw = file.path(dird, '11_qc', yid)
if(!dir.exists(dirw)) system(sprintf("mkdir -p %s", dirw))

• 设置 yid 为 'ca19a5'，这是该分析的样本标识符。
• 创建并检查工作目录路径 dirw，如果目录不存在，则使用 system(sprintf("mkdir -p %s", dirw)) 创建它。

第二部分：读取样本信息并过滤

ref = t_cfg %>% filter(yid == !!yid) %>% pull(ref)
th = rnaseq_sample_meta(yid)
tt = rnaseq_mapping_stat(yid)
res = rnaseq_cpm_raw(yid)
th = res$th; tm = res$tm; tl = res$tl; th_m = res$th_m; tm_m = res$tm_m
sum_stat_tibble(tt)

• 从 t_cfg 配置文件中根据 yid 筛选对应的参考信息。
• 使用 rnaseq_sample_meta() 函数读取样本的元数据 th。
• 使用 rnaseq_mapping_stat() 函数读取RNA-Seq映射统计信息 tt。
• 使用 rnaseq_cpm_raw() 函数读取RNA-Seq的原始CPM数据（每百万映射读数）并赋值给多个变量（如 th, tm, tl 等）。
• sum_stat_tibble(tt) 计算并显示映射统计数据的摘要。

第三部分：筛选满足条件的样本

sids_keep = tt %>% filter(mapped > 3) %>% pull(SampleID)
sum_stat_tibble(tt %>% filter(SampleID %in% sids_keep))

• 从映射统计数据 tt 中筛选出已映射的读数大于3的样本，存储其样本ID（SampleID）在 sids_keep 中。
• 输出这些筛选后样本的映射统计信息。

第四部分：修正和保存样本元数据

th2 = th %>% filter(SampleID %in% sids_keep)
th = th2
tt = tt %>% filter(SampleID %in% th$SampleID)
fh = file.path(dirw, 'meta.tsv')
write_tsv(th, fh, na='')

• 根据筛选出的样本ID，更新 th（样本元数据）。
• 根据修正后的 th，更新 tt（映射统计数据）。
• 将更新后的样本元数据保存为 meta.tsv 文件。

第五部分：计算和添加标签列

res = rnaseq_cpm(yid)
th = res$th; tm = res$tm; tl = res$tl; th_m = res$th_m; tm_m = res$tm_m
th = th %>% mutate(lab = str_c(Tissue, Genotype, sep='_'))

• 调用 rnaseq_cpm(yid) 函数重新计算CPM数据，并将其分配给 th, tm, tl 等变量。
• 在 th 数据框中，创建一个新的列 lab，将 Tissue 和 Genotype 列的值合并为一个字符串，作为标签。

第六部分：聚类分析 (Hierarchical Clustering)

tw = tm %>% select(SampleID, gid, CPM) %>% mutate(CPM=asinh(CPM)) %>% spread(SampleID, CPM)
t_exp = tm %>% group_by(gid) %>% summarise(n.exp = sum(CPM >= 1))
gids = t_exp %>% filter(n.exp >= (ncol(tw)-1) * .7) %>% pull(gid)
e = tw %>% filter(gid %in% gids) %>% select(-gid)
dim(e)

cor_opt = "pearson"
cor_opt = "spearman"
hc_opt = "ward.D"
hc_title = sprintf("dist: %s\nhclust: %s", cor_opt, hc_opt)
edist <- as.dist(1 - cor(e, method = cor_opt))
ehc <- hclust(edist, method = hc_opt)
tree = as.phylo(ehc)
lnames = ehc$labels[ehc$order]

说明：

• 数据准备：

  • tw：选取样本ID（SampleID）、基因ID（gid）和对应的CPM值，并将CPM值应用反双曲线函数（asinh()）进行转换。spread() 将数据转化为宽格式，其中每列代表一个样本。

• 过滤低表达基因：

  • t_exp：按照基因ID（gid）分组，计算每个基因表达的样本数量（n.exp），即哪些样本中该基因的CPM值大于等于1。
  • gids：选取那些在至少70%样本中有表达（n.exp >= (ncol(tw) - 1) * .7）的基因。

• 计算相关性矩阵：

  • e：从数据中筛选出表达良好的基因（gids）。移除 gid 列后，e 只包含CPM值。
  • cor_opt：定义了相关性计算的方法，可以选择 "pearson" 或 "spearman"，这两种方法分别用于计算皮尔逊相关系数和斯皮尔曼等级相关系数。
  • edist：计算基于选定相关性方法（cor_opt）的距离矩阵。as.dist(1 - cor(...)) 计算相关系数的距离（1减去相关系数）。

• 层次聚类：

  • ehc：使用 hclust() 函数进行层次聚类，ward.D 方法是最常见的聚类算法之一，用于最小化类内方差。
  • tree：通过 as.phylo() 将层次聚类结果转换为树形结构（phylogenetic tree）。

总结：

这一部分的代码进行了聚类分析，首先通过计算基因的表达情况，筛选出在大多数样本中有表达的基因。然后，使用相关性矩阵计算基因之间的相关性，并通过层次聚类方法生成基因表达的树状图。

第七部分：绘制聚类树

tp = th %>% mutate(taxa = SampleID) %>%
    select(taxa, everything())
p1 = ggtree(tree, layout = 'rectangular') +
    scale_x_continuous(expand = expand_scale(0, .2)) +
    scale_y_discrete(expand = c(.01, 0))
p1 = p1 %<+% tp + geom_tiplab(aes(label = lab, color = Tissue), size = 2.5) +
    scale_color_aaas()
fo = sprintf("%s/21.cpm.hclust.pdf", dirw)
ggsave(p1, filename = fo, width = 6, height = 6)

说明：

• 准备样本标签：

  • tp：在样本数据框 th 中新增一个 taxa 列，值为 SampleID，然后选择相关列。

• 绘制树形图：

  • 使用 ggtree() 绘制层次聚类结果的树形图，layout = 'rectangular' 设置树形图的布局为矩形。
  • 调整X轴和Y轴的显示范围和扩展。
  • 使用 geom_tiplab() 添加样本标签，其中 aes(label = lab, color = Tissue) 将样本标签设置为 lab 列，并根据 Tissue 列的值进行颜色编码。

• 保存图形：

  • 将聚类树保存为PDF文件，文件名为 21.cpm.hclust.pdf，保存路径为工作目录。

小结：

这一部分代码实现了层次聚类分析，主要通过计算基因之间的相关性来进行聚类，并通过 ggtree 绘制树形图，显示样本的聚类情况。

好的，我们继续分析最后一部分代码：tSNE分析。

第八部分：tSNE分析

require(Rtsne)
tw = tm %>% select(SampleID, gid, CPM) %>% mutate(CPM = asinh(CPM)) %>% spread(SampleID, CPM)
t_exp = tm %>% group_by(gid) %>% summarise(n.exp = sum(CPM >= 1))
gids = t_exp %>% filter(n.exp >= (ncol(tw) - 1) * .7) %>% pull(gid)
tt = tw %>% filter(gid %in% gids)
dim(tt)
tsne <- Rtsne(t(as.matrix(tt[-1])), dims = 2, verbose = T, perplexity = 1,
              pca = T, max_iter = 500)

tp = as_tibble(tsne$Y) %>%
    add_column(SampleID = colnames(tt)[-1]) %>%
    inner_join(th, by = 'SampleID')
x.max = max(tp$V1)
p_tsne = ggplot(tp, aes(x = V1, y = V2)) +
    geom_text_repel(aes(x = V1, y = V2, label = lab), size = 2.5) +
    geom_point(aes(color = Tissue, shape = Tissue), size = 2) +
    scale_x_continuous(name = 'tSNE-1') +
    scale_y_continuous(name = 'tSNE-2') +
    scale_shape_manual(values = c(0:6)) +
    scale_color_aaas() +
    otheme(legend.pos = 'top.left', legend.dir = 'v', legend.title = T,
           xtitle = T, ytitle = T,
           margin = c(.2, .2, .2, .2)) +
    theme(axis.ticks.length = unit(0, 'lines')) +
    guides(fill = F)
fp = file.path(dirw, "25.tsne.pdf")
ggsave(p_tsne, filename = fp, width = 6, height = 6)

说明：

1. 数据准备：

   • tw：与聚类分析类似，首先筛选出样本ID (SampleID)、基因ID (gid) 和对应的CPM值。然后对CPM值进行反双曲线变换（asinh()），并将数据转换为宽格式，每个样本一列。
   • t_exp：与聚类分析一样，计算每个基因的表达情况，筛选出在至少70%样本中表达的基因。
   • gids：选取那些在大多数样本中有表达的基因。
   • tt：从数据中筛选出在大多数样本中有表达的基因，保留它们的CPM值。

2. tSNE降维：

   • 使用 Rtsne 包对数据进行tSNE降维，tSNE() 函数的主要参数包括：
     • dims = 2：表示将数据降至二维空间。
     • verbose = T：在运行时显示详细信息。
     • perplexity = 1：设置困惑度（Perplexity）值，用于控制tSNE算法的表现，通常选择20-50之间的值，这里设置为1，可能是为了实验或特定数据。
     • pca = T：表示在tSNE计算前先进行PCA降维。
     • max_iter = 500：设置最大迭代次数为500次，确保算法收敛。

3. 结果合并和可视化：

   • tp：将tSNE结果（tsne$Y）转换为数据框，并将样本ID添加到结果中。同时，使用 inner_join() 将样本元数据（th）与tSNE结果按样本ID合并。

   • 绘制tSNE图：

     • ggplot() 创建tSNE的散点图，aes(x = V1, y = V2) 表示tSNE的第一维和第二维分别对应 V1 和 V2。
     • geom_text_repel() 为每个点添加标签，geom_point() 绘制样本的散点。
     • 通过 scale_color_aaas() 为不同的组织（Tissue）赋予不同的颜色，scale_shape_manual() 控制不同的形状。
     • 使用 othem() 设置图形的样式和位置，theme() 调整轴的显示方式。

4. 保存结果：

   • 将绘制的tSNE图保存为PDF文件，命名为 25.tsne.pdf，存储路径为 dirw 目录。

总结：

tSNE分析是一个常用的非线性降维方法，能够将高维数据映射到低维空间（通常是二维或三维），以便更直观地观察样本之间的关系。这里通过tSNE算法对基因表达数据进行了降维处理，并根据样本的组织类型（Tissue）对结果进行了可视化。

总结：

通过上述分析，我们逐步解析了这个R代码的不同部分。主要包括：

• 数据准备与过滤。
• 层次聚类（HClust）分析与树形图可视化。
• tSNE降维分析与散点图可视化。

每一部分都有其独立的功能，整体工作流程旨在分析和可视化基因表达数据，以便深入理解样本之间的差异性。

如果你对任何部分有进一步的疑问，或者希望继续分析其他内容，请告诉我！

原文地址：https://blog.csdn.net/2302_80012625/article/details/144816380

免责声明：本站文章内容转载自网络资源，如本站内容侵犯了原著者的合法权益，可联系本站删除。更多内容请关注自学内容网（zxcms.com）！