这个脚本的重要作用是处置惩罚和分析RNA-seq数据,进行样本的聚类分析、降维可视化,并且将结果生存下来。详细来说,它包括以下几个关键步骤:
1. 读取RNA-seq数据
- 利用rnaseq_cpm_raw()和rnaseq_cpm()两个函数分别加载原始的RNA-seq数据和颠末处置惩罚后的RNA-seq数据。
- 数据包含基因表达矩阵(tm)、样本元数据(th)以及一些基因注释信息(tl)。
2. 样本信息的加载与归并
- 从samples.xlsx文件中读取样本信息,并与RNA-seq数据中的样本元数据(th)归并。
- 将样本的基因型、处置惩罚范例、冷耐心等信息拼接成一个新的标签(lab),用于后续分析和可视化。
- 对RNA-seq数据中的表达量进行变更(利用反双曲正弦变更 asinh(CPM)),以平滑表达数据并减少大数值的影响。
3. 聚类分析(层次聚类)
- 利用层次聚类(hclust)分析样本之间的相似性,并根据干系性(皮尔逊或斯皮尔曼干系系数)计算基因表达的相似度。
- 生成聚类结果的可视化图(PDF格式),通过不同的干系性度量(如皮尔逊或斯皮尔曼)来比较样本的聚类情况。
4. 降维分析与可视化
- 通过**PCA(主成分分析)和t-SNE(t分布随机邻域嵌入)**方法将高维的基因表达数据降维到2D空间,帮助可视化数据中的样本之间的差别或聚类。
- 利用plot_pca()和plot_tsne()生成PCA和t-SNE的可视化图,并生存为PDF文件。
- 在降维图中,样本的颜色和形状可以代表不同的处置惩罚组或冷耐心状态,帮助辨认不同的样本群体。
5. 样本数据的清算与生存
- 利用complete_sample_list()函数对样本数据进行清算,去除不完整或无效的样本。
- 将终极的清算过的样本元数据(th2)生存为TSV格式的文件(01.meta.tsv),以便后续利用或备份。
6. 文件生存
- 通过ggsave()将各类分析结果(如聚类图、PCA图、t-SNE图等)生存为PDF文件,这些文件可以用于报告或进一步分析。
总结:
该脚本的焦点目标是处置惩罚RNA-seq数据并进行多种数据分析和可视化,包括:
- 数据预处置惩罚:归并样本信息、处置惩罚缺失值和进行数据变更。
- 聚类分析:通过层次聚类评估样本间的相似性。
- 降维与可视化:通过PCA和t-SNE展示样本之间的关系。
- 结果生存:将聚类和降维的可视化结果生存为文件,并将清算后的样本元数据生存为TSV文件。
这些步骤为后续的生物学分析(如基因表达模式的分析、差别分析等)奠定了基础。
第一部门:设置工作目次并加载数据
- genome = 'Osativa'
- t_cfg = read_projects(genome)
- t_cfg %>% print(n=50)
- yid = 'mo20a3'
- dirw = glue("{dird}/11_qc/{genome}/{yid}")
- if(!dir.exists(dirw)) dir.create(dirw)
复制代码 解析:
- genome = 'Osativa'
设置了基因组的名称,这里选择的是‘Osativa’(应该是指水稻基因组)。
- t_cfg = read_projects(genome)
read_projects()函数用于读取与指定基因组干系的项目配置文件,并将其存储在t_cfg中。通常这里会包含数据的元信息(例如样本、条件等)。
- t_cfg %>% print(n=50)
打印t_cfg的前50行信息,通常用于检察该基因组下全部项目标概要。
- yid = 'mo20a3'
设置实验的样本ID(yid),代表一个特定的实验或处置惩罚。
- dirw = glue("{dird}/11_qc/{genome}/{yid}")
利用glue()函数生成一个文件路径,该路径将存放QC(质量控制)结果。dird是之前定义的一个变量,大概表示根目次或数据文件夹。
- if(!dir.exists(dirw)) dir.create(dirw)
查抄是否存在上面创建的路径dirw,如果没有,则创建该目次。
这一部门的重要功能是初始化工作情况并加载数据配置文件,之后为后续的分析步骤设置目次结构。
第二部门:读取原始数据并归并样本信息
- #{{{ raw: read in
- res = rnaseq_cpm_raw(yid, genome)
- th = res$th; tm = res$tm; tl = res$tl; th_m = res$th_m; tm_m = res$tm_m
- th2 = read_xlsx(glue("~/projects/cold/data/samples.xlsx")) %>%
- replace_na(list(Name='?', ColdTolerant='?'))
- th = res$th %>% left_join(th2, by="Genotype") %>%
- mutate(lab = str_c(Genotype, Name, Treatment, ColdTolerant, sep='_'))
- tm = res$tm %>% filter(SampleID %in% th$SampleID) %>%
- mutate(value=asinh(CPM))
- #}}}
复制代码 解析:
- res = rnaseq_cpm_raw(yid, genome)
通过调用rnaseq_cpm_raw()函数,加载原始的RNA-seq数据。返回的res包含了几个重要数据框:
- th:包含样本的元数据(如样本ID、处置惩罚组等)。
- tm:包含每个样本的基因表达数据(每个基因的CPM值)。
- tl:包含基因注释信息。
- th_m、tm_m:是针对特定的处置惩罚或条件进行过计算的干系数据。
- th2 = read_xlsx(glue("~/projects/cold/data/samples.xlsx"))
从Excel文件~/projects/cold/data/samples.xlsx读取样本信息表格,th2是读取的表格数据。
- replace_na(list(Name='?', ColdTolerant='?'))
利用replace_na()函数,将缺失的样本名称(Name)和冷耐心标识(ColdTolerant)添补为问号?。
- th = res$th %>% left_join(th2, by="Genotype")
利用left_join()将res$th中的元数据和从Excel读取的th2表格数据进行归并,归并的依据是Genotype字段(基因型)。如允许以把样本的基因型信息与额外的形貌性数据联合起来。
- mutate(lab = str_c(Genotype, Name, Treatment, ColdTolerant, sep='_'))
创建一个新的列lab,将Genotype、Name、Treatment和ColdTolerant列的内容拼接在一起,并用下划线分隔。这列通常作为标签,用于标识样本。
- tm = res$tm %>% filter(SampleID %in% th$SampleID)
在表达数据res$tm中筛选出与归并后的样本信息th中SampleID匹配的行,确保只保留在th中列出的样本。
- mutate(value=asinh(CPM))
对表达矩阵中的CPM值应用反双曲正弦(asinh)变更。这是常见的对基因表达数据进行平滑处置惩罚的方式,尤其是在CPM值有较大差别时。
这一部门的关键任务是:
- 加载和归并样本元数据与表达数据。
- 对表达数据进行须要的预处置惩罚(如取对数变更)。
第三部门:聚类分析与可视化(hclust & tSNE)
- #{{{ raw: hclust & tSNE
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'pearson', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '11.hclust.p.pdf'), p1, width = 6, height = 8)
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'spearman', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '11.hclust.s.pdf'), p1, width = 6, height = 8)
- p2 = plot_pca(tm, th, pct.exp = .7, pca.center = T, pca.scale = F,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'bottom.left', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '11.pca.pdf'), p2, width = 6, height = 6)
- p3 = plot_tsne(tm, th, pct.exp = .8, perp = 3, iter = 1200, seed = 2,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'top.right', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '11.tsne.pdf'), p3, width = 6, height = 6)
- #}}}
复制代码 解析:
- 层次聚类(hclust):
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'pearson', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '11.hclust.p.pdf'), p1, width = 6, height = 8)
复制代码
- plot_hclust()函数用于进行层次聚类分析。
- tm是基因表达矩阵,th是样本信息。
- pct.exp = .7 表示只利用表达量超过70%的基因进行聚类分析。
- cor.opt = 'pearson' 选择利用皮尔逊干系系数来计算基因间的相似性。
- var.col = 'Treatment' 用样本的“Treatment”字段来为不同样本上色,以便在聚类结果中区分。
- expand.x = .4 控制图像的x轴扩展程度。
- 最后利用ggsave()将结果生存为PDF文件。
之后的部门:
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'spearman', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '11.hclust.s.pdf'), p1, width = 6, height = 8)
复制代码 这段代码执行了雷同的聚类分析,不过利用了斯皮尔曼干系系数(cor.opt = 'spearman')来代替皮尔逊干系系数。
- 主成分分析(PCA):
- p2 = plot_pca(tm, th, pct.exp = .7, pca.center = T, pca.scale = F,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'bottom.left', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '11.pca.pdf'), p2, width = 6, height = 6)
复制代码
- plot_pca()函数执行PCA分析,减少数据维度并可视化。
- pct.exp = .7 表示利用70%的基因表达数据进行分析。
- pca.center = T 表示对数据进行中心化(均值为0)。
- pca.scale = F 表示不对数据进行缩放(标准化)。
- var.shape = 'ColdTolerant' 和 var.col = 'Treatment' 用不同的形状和颜色表示不同的冷耐心和处置惩罚组。
- var.lab = 'Name' 显示样本的标签(样本名称)。
- pal.col = 'aaas' 利用指定的颜色调色板。
- 最后,利用ggsave()生存PCA图。
- t-SNE分析:
- p3 = plot_tsne(tm, th, pct.exp = .8, perp = 3, iter = 1200, seed = 2,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'top.right', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '11.tsne.pdf'), p3, width = 6, height = 6)
复制代码
- plot_tsne()用于执行t-SNE降维,并将高维数据映射到二维空间。
- pct.exp = .8 利用80%的基因表达数据。
- perp = 3 设置t-SNE的“狐疑度”参数,这影响数据点的聚集程度。
- iter = 1200 设置t-SNE迭代的次数。
- seed = 2 设置随机种子,确保结果可重复。
- var.shape = 'ColdTolerant' 和 var.col = 'Treatment' 根据样本的冷耐心和处置惩罚组着色与形状区分样本。
- legend.pos = 'top.right' 设置图例的位置。
- 最后,生存t-SNE图像。
总结:
这一部门的代码实现了:
- 层次聚类分析(利用皮尔逊和斯皮尔曼干系系数)。
- PCA降维分析,可用于观察样本在主成分空间中的分布。
- t-SNE降维分析,通过非线性方法将高维数据嵌入到二维空间中,帮助明白样本之间的相似性。
这些可视化有助于研究者了解数据的整体结构、样本的聚类情况以及不同处置惩罚组或冷耐心状态下的差别。
好的,接下来是第四部门:数据处置惩罚与生存(过滤样本信息与生存元数据)。
第四部门:数据处置惩罚与生存(过滤样本信息与生存元数据)
- #{{{ raw: filter/fix samples
- th2 = res$th
- th2 = complete_sample_list(th2)
- fh = file.path(dirw, '01.meta.tsv')
- write_tsv(th2, fh, na='')
- #}}}
复制代码 解析:
- th2 = res$th
将之前加载的样本信息 res$th 赋值给 th2。th2 存储的是包含全部样本元数据的数据框。
- th2 = complete_sample_list(th2)
complete_sample_list() 是一个自定义的函数,用于过滤或修复样本信息,通常目标是确保每个样本的元数据完整且没有缺失值。这个步骤可以用来清算样本数据,去除不完整或不符合标准的样本。
- fh = file.path(dirw, '01.meta.tsv')
创建一个文件路径 fh,用于存储样本的元数据(th2)。路径是在先前定义的工作目次 dirw 下,文件名为 01.meta.tsv。
- write_tsv(th2, fh, na='')
将处置惩罚后的 th2 数据框(样本元数据)生存为一个 TSV 文件。write_tsv() 函数将数据框写入文件,并且将缺失值(NA)生存为空字符串。
总结:
这一部门的重要任务是:
- 过滤样本数据:通过 complete_sample_list() 函数清算样本元数据,确保数据的完整性。
- 生存元数据:将清算后的样本信息生存为 TSV 格式的文件(01.meta.tsv),以便后续利用或分析。
好的,最后是第五部门:读取RNA-seq数据并重新进行聚类与可视化。
第五部门:读取RNA-seq数据并重新进行聚类与可视化
- #{{{ read in
- res = rnaseq_cpm(yid)
- th = res$th; tm = res$tm; tl = res$tl; th_m = res$th_m; tm_m = res$tm_m
- th2 = read_xlsx(glue("~/projects/cold/data/samples.xlsx")) %>%
- replace_na(list(Name='?', ColdTolerant='?'))
- th = res$th %>% left_join(th2, by="Genotype") %>%
- mutate(lab = str_c(Genotype, Name, Treatment, ColdTolerant, sep='_'))
- tm = res$tm %>% filter(SampleID %in% th$SampleID) %>%
- mutate(value=asinh(CPM))
- #}}}
- #{{{ hclust & tSNE
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'pearson', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '21.hclust.p.pdf'), p1, width = 6, height = 8)
- p1 = plot_hclust(tm, th, pct.exp = .7, cor.opt = 'spearman', var.col = 'Treatment', expand.x = .4)
- ggsave(file.path(dirw, '21.hclust.s.pdf'), p1, width = 6, height = 8)
- p2 = plot_pca(tm, th, pct.exp = .7, pca.center = T, pca.scale = F,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'bottom.left', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '21.pca.pdf'), p2, width = 6, height = 6)
- p3 = plot_tsne(tm, th, pct.exp = .8, perp = 3, iter = 1200, seed = 2,
- var.shape = 'ColdTolerant', var.col = 'Treatment', var.lab = 'Name',
- legend.pos = 'top.right', legend.dir = 'v', pal.col = 'aaas')
- ggsave(file.path(dirw, '21.tsne.pdf'), p3, width = 6, height = 6)
- #}}}
复制代码 解析:
- 重新加载RNA-seq数据:
- res = rnaseq_cpm(yid)
- th = res$th; tm = res$tm; tl = res$tl; th_m = res$th_m; tm_m = res$tm_m
复制代码 这一行代码再次加载了RNA-seq数据,执行的是 rnaseq_cpm() 函数,通常该函数会返回一系列数据框:
- th:样本的元数据(例如样本ID、处置惩罚信息等)。
- tm:基因的表达矩阵(每个基因在各个样本中的表达量)。
- tl:基因注释信息。
- th_m 和 tm_m:大概是颠末一些处置惩罚或调整的数据。
- 重新加载样本信息并归并:
- th2 = read_xlsx(glue("~/projects/cold/data/samples.xlsx")) %>%
- replace_na(list(Name='?', ColdTolerant='?'))
- th = res$th %>% left_join(th2, by="Genotype") %>%
- mutate(lab = str_c(Genotype, Name, Treatment, ColdTolerant, sep='_'))
- tm = res$tm %>% filter(SampleID %in% th$SampleID) %>%
- mutate(value=asinh(CPM))
复制代码
- 从Excel文件重新读取样本信息,确保样本的“Name”和“ColdTolerant”列中的缺失值被更换为?。
- 利用left_join()将样本信息th2与res$th按“Genotype”列归并,创建新的lab列,拼接了样本的基因型、名称、处置惩罚信息和冷耐心信息。
- 筛选出与th中样本ID匹配的基因表达数据tm,并对CPM值进行反双曲正弦变更。
- 重新进行聚类分析与可视化:
这一部门的分析和可视化步骤与前面相似,都是对基因表达数据进行聚类分析和降维可视化:
- 层次聚类(hclust):利用皮尔逊和斯皮尔曼干系系数进行层次聚类,并生存为hclust.p.pdf和hclust.s.pdf。
- PCA分析:通过plot_pca()进行主成分分析,生存为pca.pdf。
- t-SNE分析:通过plot_tsne()进行t-SNE降维,生存为tsne.pdf。
这些分析和可视化结果用于展示数据的高维结构,帮助研究人员明白不同样本之间的关系和差别。
总结:
这一部门的代码做了以下几件事情:
- 重新加载RNA-seq数据,并更新样本的元数据。
- 对基因表达数据进行处置惩罚(如变更和归并)。
- 重新进行聚类分析和降维分析,并生存可视化结果。
这与前面处置惩罚的流程雷同,只是数据加载和分析的操作略有不同。重要是为了确保数据的同等性和正确性,随后生成新的聚类与降维结果。
免责声明:如果侵犯了您的权益,请联系站长,我们会及时删除侵权内容,谢谢合作!更多信息从访问主页:qidao123.com:ToB企服之家,中国第一个企服评测及商务社交产业平台。 |