MeSH加Cytoscape你也可以绘制超高颜值富集图

MeSH加Cytoscape你也可以绘制超高颜值富集图

原创

生信技能树



生信技能树


3月19日

收录于话题
#​cytoscape的十大插件

12个

最近实验室同学在组会上分享了一篇很有意思的文章,是于  January 2021,  发表在CELL杂志的文章《Spliceosome-targeted therapies trigger an antiviral immune response in triple-negative breast cancer》,链接是:https://doi.org/10.1016/j.cell.2020.12.031

里面的GO/KEGG等生物学数据库富集注释被替换成为了 MeSH terms 的富集,而且结果还很酷炫,如下所示:

酷炫的富集分析

我看了看,文章里面的关于这个图的图例写的是:B) MeSH term enrichment analysis of top 50 resistance candidates. Enriched MeSH terms (FDR <0.1) grouped by related function. Node size represents number of shRNAs that significantly conferred resistance (R4 significant shRNAs highlighted in yellow).

也就是说,每一个基因都有一个值,就是它被多少个shRNAs所靶向,这个文章使用的是forward genetic screen with a short hairpin RNA (shRNA) library (18,370 shRNAs targeting 1,837 genes) ,它根据基因的shRNAs靶向情况把基因排序后,选取top50个基因,进行功能注释。

50个基因在附件

如下所示:

50个基因在附件

那么MeSH是何方神圣呢,我在公众号《 Y叔叔  YuLabSMU  2018-11-01》看到了这个:文章发表:Using meshes for MeSH term enrichment and semantic analyses,提示我们MeSH是个体量比GO还大的生物医学注释库,嘱咐总盯着GO看,有好多很好的注释库,都没什么人在用,比如MeSH就可以给大家提供不同的角度,挖掘不同的信息,不防试试。

步骤如下:

  • Analysis of over-represented Medical Subject Headings (MeSH) was performed using the R package ‘‘meshes’’ (v1.8.0) (Yu, 2018) with the following parameters: MeSHDb = ‘MeSH.Hsa.eg.db’, database = ‘gene2pubmed’ and category = ‘C’.
  • For Cytoscape visualization, individual genes in MeSH term categories were set as nodes and common MeSH terms as edges.

也就是说,首先你得下载和安装meshes这个包,然后使用里面的函数进行MeSH terms 的富集。

作者的MeSH富集结果也在附件

Table S2. Resistance Candidate MeSH Enrichment, Related to Figure 2

MeSH富集结果也在附件

我这里就不给出来代码了,作为一个学徒作业,反正这50个基因大家可以去文章附件复制粘贴拿到,然后R包meshes呢自己摸索一下,做一个富集。我还没有来得及看这个包,不过我猜想应该是跟Y叔的其它包其它函数用法类似,大家都是成熟的生信工程师了,该学会自己看文档了!加油!

如果是普通的富集分析

就可以使用下面的代码,当然了,也是有不同的数据库可以选择:

# 自己想办法读取 文章附件的50个基因,并且针对symbol转为ENTREZID
library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)
# top50基因
df <- bitr(unique(top50), fromType = "SYMBOL",
           toType = c( "ENTREZID" ,'GENENAME'),
           OrgDb = org.Hs.eg.db)
head(df)  
# 存储为 df 这个对象 
gene_up=df$ENTREZID
enrichKK <- enrichKEGG(gene         =  gene_up,
                       organism     = 'hsa',
                       #universe     = gene_all,
                       pvalueCutoff = 0.1,
                       qvalueCutoff =0.1)
head(enrichKK)[,1:6] 
dotplot(enrichKK)
enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
enrichKK 

比如文章的下图,就选择了 Reactome  这个数据库:

Reactome数据库注释结果

(A) Immune signaling pathway components are enriched among genes that confer resistance to spliceosome inhibition. StringDB functional enrichment analysis was performed on the top 50 SD6 resistance candidates identified by the shRNA screen. Functionally enriched Reactome Pathways with FDR < 0.001 are shown.

但是图例里面又出来了 StringDB  这个新花样!

还有一个Cytoscape的可视化

其中主要是 Table S2. Resistance Candidate MeSH Enrichment  里面的结果的解析,本身呢,Y叔的包也是提供类似的网路可视化,就是一般人想通过代码来提高颜值确实是困难,那么成熟的工具比如Cytoscape就是一个很好的选择:

至少我一般都是看了看官方文档而已:

#(3)可视化
#条带图
par(mfrow=c(2,1))
barplot(enrichKK,showCategory=20)
#气泡图
dotplot(enrichKK)
#下面的图需要映射颜色,设置和示例数据一样的geneList
# geneList = deg$logFC
# names(geneList)=deg$ENTREZID
# geneList = sort(geneList,decreasing = T)
#(3)展示top5通路的共同基因,要放大看。
#Gene-Concept Network  
cnetplot(enrichKK)

做出来的图一般来说颜值都不咋地,不过仍然是强推Y叔clusterProfiler的一些可视化方法,函数列表:

  • barplot
  • cnetplot
  • dotplot
  • emapplot
  • gseaplot
  • goplot
  • upsetplot

好几个都是以前没有介绍过的,有趣的是我准备浏览这些可视化函数的帮助文档的时候,看到了这样的话:

重点来了,Y叔特意为其包写了一本书来介绍其用法。Please go to https://yulab-smu.github.io/clusterProfiler-book/ for the full vignette.

点击数:0

美学贴面视频一例。备牙还是不备牙???-徐永钰的博客-KQ88口腔博客

美学贴面视频一例。备牙还是不备牙???

2018-11-04 15:40 阅读:2000

案例简述

     本病例的 12    13   22   23  是腭倾位的牙 ,11   21  稍稍有点前突,咬合是重度深覆盖,所以不需要考虑咬合的问题。固备牙量不用太多,正好通过贴面的厚度来补偿了 11 12   21   22之间的突度落差。患者反而觉得牙齿唇侧维度变圆润了,牙齿看上去更齐了。

     不备牙贴面适应症毕竟是少的,可以少备,但是一点不备肯定不现实,备成自切端依次向颈部逐渐变薄的样子颈部微微有大约0.2左右的肩台,就目前的材料而言其实颈部是可以最薄做到0.2毫米的。然后粘接后通过抛光这个厚度就能实现无限薄。这也是好多病例可以做的非常小,非常薄的原因。

   我们每天都在快乐的进步着,希望本视频可以带给大家一定的思路,有不足之处还望不吝赐教,谢谢

已有0人推荐

喜欢这篇文章的话,就点一下给博主推荐一下吧!

登 录

登录后才可发表评论»

最新评论

还没有评论,快来第一个评论吧!

点击数:0

单个基因多个基因的可视化

这是ggplot2可视化专题的第二个实例操作

【ggplot2的基本思路见前文总论】:基于ggplot2的RNA-seq转录组可视化:总述和分文目录

【ggplot2绘图具体策略第一篇】:测序结果概览:基因表达量rank瀑布图,高密度表达相关性散点图,注释特定基因及errorbar的表达相关性散点图绘制

【ggplot2绘图具体策略第三篇】:配对样本基因表达趋势:优化前后的散点连线图+拼图绘制

在我们获得转录组的测序结果并进行数据处理-差异分析一整套流程后,我们不仅关注整个转录组的表达趋势,还包括了特定基因在不同处理组/特征间的表达情况。这时依旧是各类ggplot2统计图的主场。

本文将要介绍的图表类型展示如下:包括分组小提琴图(B)、分组(A)/分面柱状图(C)、单基因蜂群点图拼图(D)。想必大家都知道这些图表在生物医学文献中半壁江山的地位。用R-ggplot2来绘制它们,简单明了美观,也不需要对数据进行跨平台保存读取,可谓方便。

the panels that we are going to draw in this figure

我们依旧使用通过第一、二篇介绍的整合好并差异分析过的TCGA白种人LUSC肺鳞癌mRNA-seq转录组表达数据。

数据获取

基于第一篇文章从TCGA数据库下载并整合清洗高通量肿瘤表达谱-临床性状数据,我们下载并清洗TCGA数据库中white人种的LUSC肺鳞癌mRNA-seq转录组counts数据和FPKM数据。

随后根据第二篇文章TCGA数据整合后进行DESeq2差异表达分析和基于R的多种可视化对counts数据进行了基于DESeq2的差异分析。

现在假设我们已经获得:
(样本1到344为cancer,345到386为normal)
(1)resSigAll: the subset object generated by DESeq2 containing info about all differentially expressed genes.
(2)clinical_trait: the data frame containing submitter_id and tumor_stage of selected TCGA LUSC samples. It will be used for grouping.
(3)condition_table: the data frame defining the sample_type and recording the TCGA_IDs and submitter_id of each TCGA sample. It will be used for grouping.
(4)expr_vst: vst transformed normalized counts matrix of genes and samples generated by DESeq2. It is the raw material of downstream visualization.

需要R包

ggplot2 (作图),reshape2包 (对数据框格式进行转制),ggsignif (统计学注释包),ggbeeswarm (蜂群图作图包),customLayout,和gridExtra(用于ggplot2对象拼图)。

install.packages('ggplot2')
install.packages('reshape2')
install.packages('ggsignif')
install.packages('ggbeeswarm')
install.packages('customLayout')
install.packages('gridExtra')
library('ggplot2')
library('reshape2')
library('ggsignif')
library('ggbeeswarm')
library('customLayout')
library('gridExtra')

1. 分组柱状图

假设现在我们关注一个目的基因在normal和cancer组间的表达分布,且希望看到每组内不同LUSC tumor stages间该基因的表达变化情况。(因为normal组也是来自于患者的癌旁组织,所以也可以分组)

则绘图需要:
(1)normal/cancer分组以及组内tumor stage亚组(x轴分组依据)
(2)亚组的vst counts均值(y轴,因为作柱状图,不需要每个散点的信息)
(3)亚组的vst counts标准差(errorbar)

我们随机从差异基因中选取一个基因作为可视化对象。自定义函数,输入stage i/stageii/>= stage iii可通过clinical_trait和condition_table dataframe分别匹配相应的normal和cancer样本,并分别计算均值和标准差。最后汇总得到

set.seed(10)
gene_chosen <- sample(resSigAll@rownames,1)
gene_counts_all <- expr_vst[gene_chosen,]
mean_sd_by_clinical <- function(input_stage){
  if(input_stage == '>= stage iii'){
    samples_indices <- clinical_trait$submitter_id[grep(clinical_trait$tumor_stage, pattern = 'stage iii[ab]?|stage iv')]
  } else {
    samples_indices <- clinical_trait$submitter_id[grep(clinical_trait$tumor_stage, pattern = paste0(input_stage, '[ab]?$'))]
  }
  normal_samples <- condition_table[condition_table$submitter_id %in% samples_indices & condition_table$sample_type == 'normal','TCGA_IDs']
  cancer_samples <- condition_table[condition_table$submitter_id %in% samples_indices & condition_table$sample_type == 'cancer','TCGA_IDs']
  data.frame('sample_type'=c('normal','cancer'), 'stages'=c(rep(input_stage,2)),
             'means'=c(mean(gene_counts_all[normal_samples]), mean(gene_counts_all[cancer_samples])),
             'sds'=c(sd(gene_counts_all[normal_samples]), sd(gene_counts_all[cancer_samples])))
}
mean_sd_table_by_clinical <- do.call(rbind, lapply(c('stage i', 'stage ii', '>= stage iii'), mean_sd_by_clinical))
mean_sd_table_by_clinical$sample_type <- factor(mean_sd_table_by_clinical$sample_type, 
                                                levels = c('normal','cancer'), ordered = T)

整理得到的数据框如下:

the integrated data for visualization

然后进行统计显著性框的标注。这里打算用geom_segment进行线条绘制,所以每个统一方向的线条都需要注明x/y轴的起点和终点坐标。

注:这里s指start,e指end,l指left,r指right,h指horizontal,u指upper,uh指upper horizontal。为了画出显著性注释框,这里将左竖线、右竖线、封顶横线、组间竖线分别用同一函数绘制,最后组间封顶横线用一函数绘制。

the guide for how we draw the annotation
annotation_table <- data.frame('xsl'=c(0.7, 1.7), 'xel'=c(0.7, 1.7), 'ysl'=c(7.7, 8.7), 'yel'=c(8.2, 9.2),
                               'xsr'=c(1.3, 2.3), 'xer'=c(1.3, 2.3), 'ysr'=c(7.7, 8.7), 'yer'=c(8.2, 9.2),
                               'xsh'=c(0.7, 1.7), 'xeh'=c(1.3, 2.3), 'ysh'=c(8.2, 9.2), 'yeh'=c(8.2, 9.2),
                               'xsu'=c(1, 2), 'xeu'=c(1, 2), 'ysu'=c(8.2, 9.2), 'yeu'=c(9.7, 9.7))
annotation_table2 <- data.frame('xsuh'=1 ,'xeuh'=2, 'ysuh'=9.7, 'yeuh'=9.7)

接下来就是ggplot2绘图:

◈ 将整理的mean_sd_table_by_clinical传入ggplot函数,以主要分组依据normal/cancer为x轴,y为counts均数,【fill填色依据为组内亚组tumor_stage】。

◈ geom_col函数绘制柱状图,组内亚组的柱间间距为0.7 (position),柱宽为0.5 (width),柱的边框为黑色 (color) 【fill是填充颜色,color是点/线颜色】,边框为实线 (linetype = 1)。

◈ 手动设置亚组的填充颜色。

◈ 设置误差棒的位置 (x继承ggplot分组位置,需补充指定ymax, ymin),顶端横线宽度 (width),和亚组的位置 (position: 注意与主绘图函数的间隔保持一致即可)。

◈ 不继承主参数,进行注释框的绘制 (geom_segment()),以及进行单独一行FDR值的注释 (annotate())。

◈ 设定y轴显示范围、坐标轴小标题和主题。

stage_expr <- ggplot(mean_sd_table_by_clinical, aes(x=sample_type, y=means, fill=stages))+
  geom_col(position = position_dodge(width = 0.7), width = 0.5, color='black', linetype=1)+
  scale_fill_manual(values = c('stage i'='gray20', 'stage ii'='gray50', '>= stage iii'='gray80'))+
  geom_errorbar(aes(ymin=means-sds, ymax=means+sds), width=0.1, position = position_dodge(width = 0.7), size=0.5)+
  geom_segment(inherit.aes = F, data = annotation_table, aes(x=xsl,xend=xel,y=ysl,yend=yel))+
  geom_segment(inherit.aes = F, data = annotation_table, aes(x=xsr,xend=xer,y=ysr,yend=yer))+
  geom_segment(inherit.aes = F, data = annotation_table, aes(x=xsh,xend=xeh,y=ysh,yend=yeh))+
  geom_segment(inherit.aes = F, data = annotation_table, aes(x=xsu,xend=xeu,y=ysu,yend=yeu))+
  geom_segment(inherit.aes = F, data = annotation_table2, aes(x=xsuh,xend=xeuh,y=ysuh,yend=yeuh))+
  annotate('text', label='FDR<0.001', color='black', x=1.5, y=10.1, fontface='bold', size=3)+
  ylim(0,10.3)+xlab(gene_chosen)+ylab('vst transformed normalized counts')+
  theme_bw()+
  theme(panel.grid = element_line(color = 'white'))
stage_expr

即可得到图片。注明:这里只是教学显著性注释框画法,并未分亚组分别进行假设检验。

the barplot containing sub-groups

2. 分组小提琴图

我们可能同时关注多个基因,想要同时展示它们在组间的表达情况。此时也可以类似地,将多个基因作为主要分组,而将sample_type作为fill填色依据的组内亚组,合并在一个panel进行展示。

则绘图需要:
(1)基因分组和每个基因组内的normal/cancer亚组(x轴分组和fill填色依据)
(2)每个样本相应基因的表达值(这里是vst counts)(y轴)

注:需注意不同绘图类型要求的数据形式。如作柱状图只需要每个亚组的均值,而散点图、小提琴图等需要亚组内所有样本的表达值(要进行密度展示)。又如手动绘制errorbar需提前准备每组的sd值,而绘制箱线图时不需提供分位数,ggplot2会使用输入表达数据自动计算并绘图。

这里在所有差异基因中随机选取4个基因进行展示。首先依附于expr_vst原本的数据形式整理成宽形式表格,即每个基因的表达值占据一列,并另有sample_type列进行normal/cancer分组:

format demonstration of data frame for_multigroup_vis

使用reshape2包的melt函数将数据整理成长形式,此时主分组和亚组信息就可以直接分别从genes和sample_type列传递给ggplot了。

format demonstration of data frame melted_var_table

然后将主分组genes进行factor化排序。

set.seed(29)
for_multigroup_vis <- t(expr_vst[sample(resSigAll@rownames, 4),])
for_multigroup_vis <- data.frame(for_multigroup_vis, 'sample_type'=condition_table$sample_type, stringsAsFactors = F)
melted_var_table <- reshape2::melt(for_multigroup_vis, id.vars='sample_type', 
                                   variable.name='genes', value.name='vst_counts')
#tell ggplot2 the order of graphing within each gene group.
melted_var_table$sample_type <- factor(melted_var_table$sample_type, levels = c('normal','cancer'), ordered = T)

和前面的示例相同,建立坐标表以使用geom_segment函数绘制各组的统计显著性框,并建立坐标和标签表以使用geom_text函数添加各组FDR文字注释。

FDR_text_table <- data.frame('gene'=as.character(unique(melted_var_table$genes)),'start'=c(0.8,1.8,2.8,3.8), 'end'=c(1.19, 2.19,3.19,4.19),
                                'height'=c(22,15,16.5,16), annotations=c('FDR=2.15e-91','FDR=3.40e-59','FDR=1.74e-89','FDR=1.07e-48'))
segment_table <- data.frame('xsl'=c(0.8,1.8,2.8,3.8), 'xel'=c(0.8,1.8,2.8,3.8), 'ysl'=c(17,10,15.3,14.7), 'yel'=c(21,14,15.5,15),
                            'xsr'=c(1.2,2.2,3.2,4.2), 'xer'=c(1.2,2.2,3.2,4.2), 'ysr'=c(20.5,13.8,14,14.5), 'yer'=c(21,14,15.5,15),
                            'xsh'=c(0.8,1.8,2.8,3.8), 'xeh'=c(1.2,2.2,3.2,4.2), 'ysh'=c(21,14,15.5,15), 'yeh'=c(21,14,15.5,15))

ggplot2绘图的思路基本相同,不再详细描述。值得注意的点是:

◈ 此次映射数据时以genes为x轴主分组,sample_type为fill填色亚分组,每个样本的vst_counts映射至y轴。

◈ 使用geom_violin()绘制小提琴图,position参数设置主组内亚组间距离,使用scale设置组间面积恒定 ‘area’。(或设置为width/count,分别表示每组有相同的最大小提琴宽度和宽度与观察值数成比例)

#draw the plot.
multiple_vio <- ggplot(melted_var_table, aes(x=genes, y=vst_counts, fill=sample_type))+
  geom_violin(position = position_dodge(width=0.8), scale = 'area')+
  geom_boxplot(position = position_dodge(width = 0.8), 
               outlier.size = 0, width = 0.1, show.legend = FALSE)+
  scale_fill_manual(values = c('cancer'='orange','normal'='blue'))+
  geom_text(inherit.aes = F,data = FDR_text_table, aes(x=gene, y=height, label=annotations), size=2.5, fontface='bold')+
  geom_segment(inherit.aes = F,data=segment_table,aes(x=xsl,xend=xel,y=ysl,yend=yel))+
  geom_segment(inherit.aes = F,data=segment_table,aes(x=xsr,xend=xer,y=ysr,yend=yer))+
  geom_segment(inherit.aes = F,data=segment_table,aes(x=xsh,xend=xeh,y=ysh,yend=yeh))+
  theme_bw()+theme(panel.grid = element_line(colour = 'white'), plot.title = element_text(hjust = 0.5))+
  ylab('normalized_counts')+xlab('genes')+ylim(3,24)+
  ggtitle('violin plot of multiple genes')
multiple_vio

获得图形为:

the violin plot containing multiple groups and sub-groups

3. 分面柱状图

上述两个图片的策略均是二次分组,x映射主分组,fill映射亚组。另一种策略是使用ggplot2的facet.grid()函数进行分面展示,即每个主分组获得一个panel分别可视化。

则绘图需要:
(1)同时含有基因和normal/cancer信息的独立分组(x轴分组依据)
(2)每个组隶属的基因信息(x轴分面)
(3)每个组隶属的样本类型(fill填充依据)
(4)每个组的表达值均数(y轴)
(5)每个组的表达值标准差(errorbar)

整理数据,沿用之前随机选取的4个基因并整理成长格式的表达data frame。自建函数,向量化批量计算每个基因-每种样本类型独立组的表达平均值和标准差。以”gene_sample_type“的格式标记为8个独立组,并添加annotation列便于geom_text的注释。最后将独立分组因子化用于手动排序,人工将同一个基因的normal/cancer组相邻放置。

#draw the graphs by facet.
melted_var_table2 <- melted_var_table
#calculate the mean and sd of each group.
mean_sd <- lapply(as.character(unique(melted_var_table2$genes)), function(x){
  mean_normal <- mean(melted_var_table2$vst_counts[melted_var_table2$genes == x & melted_var_table2$sample_type == 'normal'])
  mean_cancer <- mean(melted_var_table2$vst_counts[melted_var_table2$genes == x & melted_var_table2$sample_type == 'cancer'])
  sd_normal <- sd(melted_var_table2$vst_counts[melted_var_table2$genes == x & melted_var_table2$sample_type == 'normal'])
  sd_cancer <- sd(melted_var_table2$vst_counts[melted_var_table2$genes == x & melted_var_table2$sample_type == 'cancer'])
  data.frame('gene'= c(x,x), 'sample_type'=c('normal','cancer'),
             'mean'=c(mean_normal, mean_cancer), 'sd'=c(sd_normal, sd_cancer))})
mean_sd <- do.call(rbind, mean_sd)
mean_sd_annotation <- data.frame(mean_sd, 'groups'=paste(mean_sd$gene, mean_sd$sample_type, sep = '_'),
                                 'annotation'= rep(c('','****')), stringsAsFactors = F)
mean_sd_annotation$groups <- factor(mean_sd_annotation$groups, levels = unique(mean_sd_annotation$groups),
                                   ordered = T)

整理好的数据框为如下,含gene分类,sample_type分类,将两列字符串粘贴在一起而形成的独立组信息,每组的mean/sd值,以及每个组对应的annotation信息:

format demonstration of the data frame mean_sd_annotation

ggplot2绘图:值得注意的点是:

◈ ggplot函数x轴映射groups,即各个独立组。fill以sample_type分别配色。

◈ geom_text补充文字的y轴位置,即y=mean+sd+0.5。

◈ 使用facet_grid()函数,以输入数据框的gene列按行分面 (~gene),每个分面只保留有数据的x轴部分 (scale=’free’)。

p <- ggplot(mean_sd_annotation, aes(x=groups, y=mean, fill=sample_type))+
  geom_col(width=0.4, linetype=1, color='black')+
  geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd), width=0.1)+
  geom_text(aes(y=mean+sd+0.5,label=annotation), size = 3, fontface='bold')+
  scale_fill_manual(values=c('normal'='blue', 'cancer'='orange'))+
  facet_grid(~gene, scales='free')+
  scale_x_discrete(label=c('normal','cancer'))+
  ylab('vst-normalized counts')+xlab('groups')+
  theme_bw()+theme(panel.grid = element_line(colour = 'white'), legend.position = 'none')
p

获得图形如下:

每个框都属于一个facet的分面柱状图

注:使用facet.grid函数时需要注意各个函数间需要有统一的映射数据。比如:本例如果使用geom_violin(),则需输入所有样本的表达值而不是均值,这会使得文字注释的变量长度与输入数据框长度不一而无法合并到一个数据框中共同传入ggplot()。若geom_text不继承全局参数,会导致facet.grid认为每个分面在x轴的每个分组上均有数据,导致布局混乱。

总之,分facet这个方法是比较局限的。

4. 拼接多个单基因蜂群图

如果我想使用更多样性的可视化方法呢?比如在样本例数较多时使用蜂群图而不是较为死板的柱状图?

an example of beeswarm plot downloaded from an article of Plos One

R中有beeswarm包,加载后使用geom_beeswarm()函数可实现本功能。但由于未知的原因,它和geom_jitter()都不能像geom_violin()一样进行二次分组,如果采取第二种作图策略直接设置ggplot(data, aes(x=gene, y=vst_counts, color=sample_type)),他们仅表现为geom_point即可完成的单纯散点图。这不是我们想要的效果。如下图:

an erroneous example

所以这里我使用产生多个单基因一次分组的蜂群散点图,并用拼图包将所有单图拼接。

绘图需要:
(1)各基因的sample_type (x轴);
(2)各组内样本的表达值 (这里是vst-counts) (y轴);
(3)各组的平均值和标准差 (自制mean±sd errorbar)。

首先整理数据框,获得各样本的gene,sample_type,vst_counts数据框,各组的均值、标准差数据框,和各基因组间显著性标记数据框。

#draw the bee plot of the same data by another manner (by joining graphs together).
melted_var_table3 <- melted_var_table
#we should manually tell the ggplot2 how we oder the x axis by transforming groups into factors.
melted_var_table3$sample_type <- factor(melted_var_table3$sample_type, 
                                        levels = c('normal','cancer'), ordered = T)
#calculate the mean and sd of each group.
mean_sd <- lapply(as.character(unique(melted_var_table3$genes)), function(x){
  mean_normal <- mean(melted_var_table3$vst_counts[melted_var_table3$genes == x & melted_var_table3$sample_type == 'normal'])
  mean_cancer <- mean(melted_var_table3$vst_counts[melted_var_table3$genes == x & melted_var_table3$sample_type == 'cancer'])
  sd_normal <- sd(melted_var_table3$vst_counts[melted_var_table3$genes == x & melted_var_table3$sample_type == 'normal'])
  sd_cancer <- sd(melted_var_table3$vst_counts[melted_var_table3$genes == x & melted_var_table3$sample_type == 'cancer'])
  data.frame('gene'= c(x,x), 'sample_type'=c('normal','cancer'),
             'mean'=c(mean_normal, mean_cancer), 'sd'=c(sd_normal, sd_cancer))})
mean_sd <- do.call(rbind, mean_sd)
#we have gotten the FDRs from the object of DE analysis before.
annotation_FDR <- data.frame('gene'=as.character(unique(melted_var_table$genes)),
                             'annotations'=c('FDR=2.15e-91','FDR=3.40e-59','FDR=1.74e-89','FDR=1.07e-48'))

ggplot2作图。由于拼图R包customLayout拼图的输入对象是多个ggplot2绘图对象的list,因此这里创建一个自定义ggplot2绘图function,之后配合lapply函数建立多个单基因蜂群图的list。

几个注意点:
◈ geom_beeswarm()绘制蜂群点图,cex设置点之间的排列距离,priority设置一个组中点的排列方式。

◈ 使用geom_point(),以均值标准差数据框作为输入,标记每个组的平均值,设置shape为3即加号 (+),方便与errorbar融合在一起。

◈ 使用geom_errorbar(),以均值标准差数据框作为输入,绘制sd值errorbar,在图中与均值的+号融合形成三分errorbar形状。

◈ 使用一个新的R包ggsignif进行统计学注释。geom_signif()仍继承全局参数,comparisons输入所有进行比较的组间名称的list (list(c(‘group1’, ‘group2’), c(‘group3’, ‘group4’), …)) ,annotations参数则是按comparisons参数中输入比较组的顺序添加文字注释。本包只能进行一次分组的比较,而不能进行由fill/color划分的亚组间的比较注释。

◈ 为了防止不同基因间y轴区间差异过大,使用ylim()设置每个图的y轴作图区间仅分布在本基因样本最小表达值-1.5,到最大表达值+1.5之间。

multi_bee_plot <- function(gene_for_vis){
  library(ggplot2)
  library(ggbeeswarm)
  library(ggsignif)
  ggplot(melted_var_table3[melted_var_table3$genes == gene_for_vis,], aes(x=sample_type, y=vst_counts, color=sample_type))+
  geom_beeswarm(cex = 1.5, size=0.5, priority = 'ascending')+ #priority=none/random/ascending/descending
  geom_point(inherit.aes = F, data=mean_sd[mean_sd$gene == gene_for_vis,],
             aes(x=sample_type, y=mean), shape=3, color='black', size=5)+
  geom_errorbar(inherit.aes = F, data=mean_sd[mean_sd$gene == gene_for_vis,],
                aes(x=sample_type, ymin=mean-sd, ymax=mean+sd), width=0.5, size=0.5)+
  scale_x_discrete(labels=c('normal','cancer'))+
  scale_color_manual(values=c('normal'='blue','cancer'='orange'))+
  geom_signif(comparisons = list(c('normal','cancer')), annotations = annotation_FDR$annotations[annotation_FDR$gene == gene_for_vis],
              color='black', size = 0.5, textsize=3, fontface='bold')+
  xlab(gene_for_vis)+ylab('vst normalized counts')+
  ylim(min(melted_var_table3[melted_var_table3$genes == gene_for_vis,'vst_counts'])-1.5,
       max(melted_var_table3[melted_var_table3$genes == gene_for_vis,'vst_counts'])+1.5)+
  theme_bw()+
  theme(panel.grid = element_line(color = 'white'), legend.position = 'none')
}

使用拼图R包customLayout进行拼图。

思路是用lay_new函数创建新画布,建立matrix,用其中的行、列、和数字分别代表画布中容纳图片的行、列数和图片排列的顺序。

lay <- lay_new(mat=matrix(1:4, ncol = 4), widths=c(2,2,2,2),heights = 1)
lay_show(lay)

lay_show函数用于展示创建的画布:

the canvas we generated

lapply函数创建绘图对象list,并用lay_grid函数进行拼图输出。

multiple_beeswarm <- lapply(as.character(unique(melted_var_table2$genes)), multi_bee_plot)
lay_grid(grobs = multiple_beeswarm, lay = lay)

获得如下组图:

beeswarm plots collage

一共四种解决方案,希望本文能对大家有所帮助!

喜欢的朋友就学起来,点个赞吧!

点击数:0

KEGG分析工具:pathview基本操作 – 简书

KEGG分析工具:pathview基本操作

32018.07.10 20:48:46字数 1,264阅读 12,857
什么是KEGG

KEGG(Kyoto Encyclopedia of Genes and Genomes)是系统分析基因功能、基因组信息数据库,它有助于研究者把基因及表达信息作为一个整体网络进行研究。KEGG将基因组信息和高一级的功能信息有机地结合起来,通过对细胞内已知生物学过程的计算机化处理和将现有的基因功能解释标准化,对基因的功能进行系统化的分析。KEGG的另一个任务是一个将基因组中的一系列基因用一个细胞内的分子相互作用的网络连接起来的过程,如一个通路或是一个复合物,通过它们来展现更高一级的生物学功能。KEGG主要包含以下数据库:

解释基因与功能之间的关系
为什么要用KEGG的代谢通路

KEGG提供的整合代谢途径(pathway)查询十分出色,包括碳水化合物、核苷、氨基酸等的代谢及有机物的生物降解,不仅提供了所有可能的代谢途径,而且对催化各步反应的酶进行了全面的注解,包含有氨基酸序列、PDB库的链接等等。KEGG是进行生物体内代谢分析、代谢网络研究的强有力工具。与其他数据库相比,KEGG 的一个显著特点就是具有强大的图形功能,它利用图形而不是繁缛的文字来介绍众多的代谢途径以及各途径之间的关系。

  • 从功能出发,研究功能到通路到基因,迅速锁定某些功能的基因;
  • 从基因出发,获得某个基因在信号通路中的角色(上下游关系)和生物学功能;
  • 发现涉及通路的差异变化和功能分布
  • 形象的图形使我们直观地对某一个基因有了一个由点及面的印象。
KEGG代谢通路怎么看

在KEGG中有两种代谢图

  • 参考代谢通路图reference pathway,是根据已有的知识绘制的概括的、详尽的具有一般参考意义的代谢图,这种图上就不会有绿色的小框,而都是无色的,所有的框都可以点击查看更详细的信息;

  • 特定物种的代谢图species-specific pathway,会用绿色来标出这个物种特有的基因或酶,只有这些绿色的框点击以后才会给出更详细的信息。

这两种图很好区分,reference pathway 在KEGG中的名字是以map开头的,比如map00010,就是糖酵解途径的参考图;而特定物种的代谢通路图开头三个字符不是map而是种属英文单词的缩写(应该就是一个属的首字母+2个种的首字母)比如酵母的糖酵解通路图,就是sce00010,大肠杆菌的糖酵解通路图就应该是eco00010吧。

代谢通路中各种符号标识 :

  • K+num:基因ID号,表示在所有同源物种中具有相似结构或功能的一类同源蛋白
  • ko+num: 代谢通路名称,表示一个特定的生物路径
  • M+ num: 模块名称
  • C+ num: 化合物名称
  • E-,-,-,-: 酶名称
  • R + num : 反应名
  • RC+ num : 反映类型
  • RP+num: 反应物对

图例作用关系:

如何绘制代谢通路图

控制通路图将自己的基因绘制在一个通路图中

  • 线分析工具—ipath2.0
  • R包:pathview
    今天简单介绍一下这个包。
#install 
source("http://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("pathview")
library("pathview")
setwd("C:\Users\Administrator\Desktop\pathview")
rm(list=ls())
#载入数据文件
data(gse16873.d)
data(demo.paths)
cpd.data <- read.csv("sim.cpd.data2.csv",row.names=1)
?pathview
View(gse16873.d)

基因表达变化数据框gse16873.d如下所示,行是基因ID,列是样本ID,变化范围是-1到1.`

代谢通路数据demo.paths结构如下:

基本参数:

#Usage
pathview(gene.data = NULL, cpd.data = NULL, pathway.id,
species = "hsa", kegg.dir = ".", cpd.idtype = "kegg", gene.idtype =
"entrez", gene.annotpkg = NULL, min.nnodes = 3, kegg.native = TRUE,
map.null = TRUE, expand.node = FALSE, split.group = FALSE, map.symbol =
TRUE, map.cpdname = TRUE, node.sum = "sum", discrete=list(gene=FALSE,
cpd=FALSE), limit = list(gene = 1, cpd = 1), bins = list(gene = 10, cpd
= 10), both.dirs = list(gene = T, cpd = T), trans.fun = list(gene =
NULL, cpd = NULL), low = list(gene = "green", cpd = "blue"), mid =
list(gene = "gray", cpd = "gray"), high = list(gene = "red", cpd =
"yellow"), na.col = "transparent", ...)

Note that gene.data and cpd.data can't be NULL simultaneously.

gene.data是需要提供的基因向量,默认是Entrez_ID。其由gene.idtype决定
cpd.data 指的药物分子的名称向量。
Pathway.id指的是在KEGG中的ID。
kegg.native默认是TRUE输出完整pathway的png格式文件,反之输出仅是输入的基因列表的pdf文件。
Map.null默认是TRUE,当使用FALSE时其pdf的文件图像会更漂亮
Split.group 主要是在kegg.native为FALSE的时候会起到一定的作用,主要是将在同一个反应的基因归在一起。
new.signature=FALSE将会将标签去掉,只显示图像

#原始的kegg.native=TRUE时的图像绘制
pv.out <- pathview(gene.data = gse16873.d[, 1], pathway.id =
                     demo.paths$sel.paths[3], species ="hsa", out.suffix = "gse168731", gene.idtype =
                     "entrez", gene.annotpkg = NULL,min.nnodes = 3, kegg.native =TRUE,
                   map.null = FALSE, expand.node =FALSE,split.group =FALSE, map.symbol =
                     TRUE,new.signature=FALSE)

#kegg.native =FALSE:
pv.out <- pathview(gene.data =gse16873.d[, 1], pathway.id =
                     demo.paths$sel.paths[3], species ="hsa", out.suffix = "gse168732", gene.idtype =
                     "entrez", gene.annotpkg = NULL,min.nnodes = 3, kegg.native =FALSE,
                   map.null = FALSE, expand.node =FALSE,split.group =FALSE, map.symbol =
                     TRUE,new.signature=FALSE)
#进一步如果想将所有同一个反应的基因归在一起,那么需要设置参数split.group =TRUE:
pv.out <- pathview(gene.data =gse16873.d[, 1], pathway.id =
                     demo.paths$sel.paths[1], species ="hsa", out.suffix = "gse168733", gene.idtype =
                     "entrez", gene.annotpkg = NULL,min.nnodes = 3, kegg.native =FALSE,
                   map.null = FALSE, expand.node =FALSE,split.group =TRUE, map.symbol =
                     TRUE,new.signature=FALSE)
#KEGG view: both gene and compound data
sim.cpd.data=sim.mol.data(mol.type="cpd", nmol=3000)

pv.out <- pathview(gene.data = gse16873.d[, 1], cpd.data = sim.cpd.data,
                   pathway.id = demo.paths$sel.paths[3], species = "hsa", out.suffix =
                     "gse16874.cpd", keys.align = "y", kegg.native = TRUE, key.pos = demo.paths$kpos1[3])
#multiple states in one graph

pv.out <- pathview(gene.data = gse16873.d[, 1:6], 
                   cpd.data = cpd.data[, 1:6], pathway.id = demo.paths$sel.paths[3], 
                   species = "hsa", out.suffix = "gse16875", keys.align = "y", 
                   kegg.native = TRUE, match.data = FALSE, multi.state = TRUE, same.layer = TRUE)

hsa00640.gse16875.multi.png

pathview在线版
Pathview: An R package for pathway based data integration and visualization
数据分析-【KEGG相关包-clusterProfiler,Pathview的学习】
经典信号通路作图工具包
R语言实现KEGG通路富集可视化
KEGG信号通路的展示
KEGG简介、如何使用KEGG进行通路富集?

点击数:0

可视化kegg通路-pathview包 | KeepNotes blog

可视化kegg通路-pathview包

Posted on

2018-08-22

Edited on
2019-11-09

In

Bioinformatics-Notes

,

Basic

Symbols count in article:
4.7k

Reading time ≈
4 mins.

Pathview是一个用于整合表达谱数据并用于可视化kegg通路的一个R包,其会先下载kegg官网上的通路图,然后整合输入数据对通路图进行再次渲染(加工?),从而对kegg通路图进行一定程度上的个性化处理 Pathview是一个bioconductor包,正常安装即可

    source("https://bioconductor.org/biocLite.R")
    biocLite("pathview")

以其自带测试数据集,对Pathview包的可视化功能做个整理,用法很简单,主要为围绕着pathview这个作图函数(参数很多,用?pathview查看下各个参数的说明),只是需要将数据整理或者整合成其需要的输入数据格式即可

先看下demo数据是什么格式,列名是每个样本名,行名是每个gene的entrez id

> head(gse16873.d)
               DCIS_1      DCIS_2       DCIS_3      DCIS_4       DCIS_5      DCIS_6
10000     -0.30764480 -0.14722769 -0.023784808 -0.07056193 -0.001323087 -0.15026813
10001      0.41586805 -0.33477259 -0.513136907 -0.16653712  0.111122223  0.13400734
10002      0.19854925  0.03789588  0.341865341 -0.08527420  0.767559264  0.15828609
10003     -0.23155297 -0.09659311 -0.104727283 -0.04801404 -0.208056443  0.03344448
100048912 -0.04490724 -0.05203146  0.036390376  0.04807823  0.027205816  0.05444739
10004     -0.08756237 -0.05027725  0.001821133  0.03023835  0.008034394 -0.06860749

先看下Pathview最常见的一种用法:将某个样本的表达量(测试数据已经是归一化后的表达量);其实也可以任何列数据,不仅仅是表达量数据,也可以是foldchange等数据,人为特定的数值型数据也行;最后以color bar的形式在kegg通路图上的对应节点展示;如下例子所示,取第一个样本的数值最为gene.data,通路选择04110,物种为hsa

p <- pathview(gene.data = gse16873.d[, 1], pathway.id = "04110", species = "hsa", out.suffix = "gse16873")
hsa04110.pathview

图中右上角color bar是可以通过参数limit来调整的,默认是limit=list(gene=1, cpd=1),即gene.data的限定范围在(-1,1),cpd.data的限定范围在c(-1,1);如果有其他需要可以设置为limit=list(gene=2, cpd=1)等;如果想最大值和最小值不对称则limit = list(gene=c(-1,2))即可

如果想将color bar分割的更加密集一些,则可以修改默认参数bins = list(gene = 10, cpd = 10)为20。。。只能说Pathview作者好细心。。

接下来看看Pathview如何展示整合数据的,如基因和代谢组的数据在kegg通路上整合可视化

先载入代谢物相关数据,如

sim.cpd.data=sim.mol.data(mol.type="cpd", nmol=3000)
> head(sim.cpd.data)
     C02787      C08521      C01043      C11496      C07111      C00031 
-1.15259948  0.46416071  0.72893247  0.41061745 -1.46114720 -0.01890809 

然后通过cpd.data参数将代谢物数据加入到pathview函数中,如:

p <- pathview(gene.data = gse16873.d[, 1], cpd.data = sim.cpd.data, 
               pathway.id = "00640", species = "hsa", out.suffix = "gse16873.cpd", 
               keys.align = "y", kegg.native = T, key.pos = "topright")
hsa00640.cpd.pathview

keys.align = "y"控制color bar对齐方式,key.pos = "topright"则是对应color bar在图上的位置(主要有bottomleft,bottomright,topleft以及topright)

Pathview不仅能整合不同组学数据(上述例子是转录组和代谢组),还可以整合多个样本的数据,比如我们在上面看到gene.data只导入了一列数据,其实Pathview支持多列(也就是多个样本)数据的导入,cpd.data也一样

如按照说明文档中的模拟一组多样本的代谢物表达谱数据(也似乎是归一化后的)

sim.cpd.data2 = matrix(sample(sim.cpd.data, 18000, replace = T), ncol = 6)
rownames(sim.cpd.data2) = names(sim.cpd.data)
colnames(sim.cpd.data2) = paste("exp", 1:6, sep = "")
> head(sim.cpd.data2, 3)
             exp1      exp2       exp3       exp4       exp5       exp6
C02787 -0.5425224 1.7940544 -0.2629972  0.2729004 -0.4897083 1.05131740
C08521 -1.1903358 0.4448658  2.6074747 -0.9163451  0.1239377 0.57827710
C01043  0.3391817 1.6855815  1.0203767 -1.3184792  0.4727454 0.03381888

然后还是使用pathview函数作图,相比上面例子增加了参数match.data = F,主要用于表示gene和cpd样本数是否要匹配;multi.state = T则表示多个样本在同一个图上显示

p <- pathview(gene.data = gse16873.d[, 1:3], cpd.data = sim.cpd.data2[, 1:2], 
               pathway.id = "00640", species = "hsa", out.suffix = "gse16873.cpd.3-2s", 
               keys.align = "y", kegg.native = T, match.data = F, multi.state = T, same.layer = T)
hsa00640.cpd.3-2s.multi.pathview

正如之前说的,Pathview不仅能展示连续性数据,还可以展示离散型数据(如:差异显著OR不显著,上调OR下调等),我先模拟一个(0,1)离散型数据

testdata2 <- data.frame(sample1=sample(c(0,1), 10000, replace = T), stringsAsFactors = F)
row.names(testdata2) <- row.names(gse16873.d)[1:10000]

然后在pathview函数中将limit = list(gene = c(0,1))bins = list(gene = 1)设置为一致;由于默认颜色是低(绿色),中(灰色),高(红色),而我们测试数据是0,1两个数字,所以需要将mid设置为红色mid = list(gene = "red"),这样0就是绿色,1就是红色了

p <- pathview(gene.data = testdata, pathway.id = "00640", species = "hsa", 
               out.suffix = "test.genes",discrete = list(gene = T), kegg.native = T, 
               limit = list(gene = c(0,1)), bins = list(gene = 1), mid = list(gene = "red"))
hsa00640.test.genes

从我们之前的测试数据中可看到,pathview函数默认使用的Id是entrez id,从默认参数gene.idtype="entrez"可得;对于绝大多数真核生物来说,其entrez id确实等于其对应的kegg id,但是有些物种还是特例的,比如拟南芥。。。当然Pathview包给出了一些id转化方法(但我一般不用,因为我比较喜欢自己来转化);因此一般我们能拿到的就已经是kegg id了,比如拟南芥的ath:AT4G17260,这是需要设置gene.idtype=kegg,这里需要注意的,gene.data参数输入的不是ath:AT4G17260而是AT4G17260,下面以一个例子说明

先模拟一个数据

testid <- data.frame(expr = rnorm(5), stringsAsFactors = F)
row.names(testid) <- c("AT4G17260", "AT2G14170", "AT2G20420", "AT1G06550", "AT1G36160")

然后还是以通路00640为例,物种是拟南芥

p <- pathview(gene.data = testid, pathway.id = "00640", species = "ath", 
               gene.idtype = "kegg", out.suffix = "testid")
ath00640.testid

但是当物种是kegg未收录的物种,这是可能用的不再是kegg id还是K号,如K00016,这时的通路则是KEGG ortholog pathways,物种设定为ko,即:species = "ko"

testko <- data.frame(expr = rnorm(5), stringsAsFactors = F)
row.names(testko) <- c("K00016", "K01026", "K00140", "K18371", "K17741")
ko00640.testko

上述大概是Pathview包主要的用法了,如果有需要可以再从头看下其说明文档http://bioconductor.org/packages/release/bioc/vignettes/pathview/inst/doc/pathview.pdf

本文出自于http://www.bioinfo-scrounger.com转载请注明出处

点击数:0

用limma包进行多组差异表达分析_今天也是个妖精头子呀的博客-CSDN博客_limma 差异表达分析

用limma包进行多组差异表达分析

置顶
今天也是个妖精头子呀
2020-03-15 00:20:16

4444

分类专栏:
生物信息学基础知识

写在前面:最近在使用limma包进行差异表达分析,参考了网上许多教程都觉得说的云里雾里,很不清楚。经过我自己一段时间非常痛苦的钻研,弄明白了,解决了我的实际需求。于是决定将我的分析经验写下来,分享给需要的人。

 首先加载前期预处理好的表达矩阵。(我的原始表达矩阵文件已附在文后,大家有需要可以下载实践)

表达矩阵格式如下图所示,行名为基因名,列名为样本名。

创建样本分类信息表。

所谓样本分类信息表,实际上就是值你的样本名一一对应的属性是什么(譬如case还是control)。

我自己根据GEO数据库中该芯片的分类信息先自己在excel中整理了一张表。因为这张芯片原作者在上传数据的时候,样本分类数据很乱,实验组和对照组的样本有时候交错分布,譬如ID为GSM2612287的样本是第192个样本,是对照组即健康人群。而ID为GSM2612128的样本是第33个样本,取样自重症抑郁症患者。我们在对样本进行差异表达分析之前,一定要做的工作是弄明白每一个样本的属性,我做这个工作,实际上是对样本的属性进行了手工注释。

然后我将上图中的第二第三列,复制单独保存在一个txt文档中。(即gene_list.txt)

第一列(1,2,3,……)只是为了辅助排序的东西,因为这个1,2,3,就是样本在总的样本中的位置,这个表格就是告诉我们第几个样本,它是对照组,第几个样本,取样自重症抑郁症患者。

1 HC
2 HC
3 HC
4 HC
…… ……

加载上述分类表gene_list.txt。并按照数据框的第一列升序排序。

取第二列,这样处理得到的就是每一个位置所对应的样本属性信息。我们最终想要得到的也就是向量group。(实际情况大家可能不会像我这样麻烦,具体问题具体分析,只要得到一个类如group的向量)

注:其中HC代表的是健康对照组,MDD指重症抑郁症;anxiety,即焦虑症。

接着,加载limma包,构造如下矩阵。(这一段可以照抄。根据自己的实际情况替换变量)

然后,我们规定哪一组数据与哪一组数据比较。

这里也是原来我比较困惑的地方,因为我的数据是有MDD,anxiety两个实验组,HC对照组。而我实际只需要比较MDD与HC的差异表达情况,而一般的教程中并未涉及或者介绍模糊。

下面这段的代码的意思就是比较design矩阵中HC与MDD level。(关于多组分析的其它疑惑,可以参考知乎文章

contrast.matrix<-makeContrasts("HC-MDD",levels=design)

完成以上两个表格之后,就可以直接跑下面的代码,得到差异表达数据了。(有时间具体分析什么意思)

分析得到的nrDEG如下图所示:

 保存分析好的差异表达数据。

write.table(nrDEG,"limma_GSE98793.txt")

 全部代码:

 熬夜写到这儿,一定会有些疏漏和错误,欢迎指正,共同进步!

【数据链接】https://me.csdn.net/download/weixin_40640700

重新修改了上传数据所需的积分,设置为零。可是系统还是自动调为了1。

如果参考数据不能顺利下载的话,可以邮箱2456392738@qq.com,我邮箱发送。

点击数:0

使用htseq-count进行定量分析 – 简书

使用htseq-count进行定量分析

0.3412020.03.10 00:11:53字数 2,314阅读 2,972

和featurecounts一样,htseq-count也是一款进行raw count定量的软件。该软件采用python语言进行开发,集成在HTseq这个包中。

对于python的包,通过pip可以方便的进行安装,代码如下

pip install HTSeq

HTSeq提供了许多处理NGS数据的功能,htseq-count只是其中进行定量分析的一个模块。

HTSeq使用注意事项

1、HTSeq是对有参考基因组的转录组测序数据进行表达量分析的,其输入文件必须有SAM和GTF文件。

2、一般情况下HTSeq得到的Counts结果会用于下一步不同样品间的基因表达量差异分析,而不是一个样品内部基因的表达量比较。因此,HTSeq设置了-a参数的默认值10,来忽略掉比对到多个位置的reads信息,其结果有利于后续的差异分析。

3、输入的GTF文件中不能包含可变剪接信息,否则HTSeq会认为每个可变剪接都是单独的基因,导致能比对到多个可变剪接转录本上的reads的计算结果是ambiguous,从而不能计算到基因的count中。即使设置-i参数的值为transcript_id,其结果一样是不准确的,只是得到transcripts的表达量。

htseq-count的设计思想和featurecounts非常类似,也包含了feature和meta-features两个概念。对于转录组数据而言,feature指的是exon, 而meta-feature可以是gene, 也可以是transcript。

进行定量分析需要以下两个文件

1、比对的BAM/SAM文件

2、基因组的GTF文件

对于双端数据,要求输入sort之后的BAM文件。

由于序列读长的限制和基因组的同源性,一条reads可能比对到多个基因上,而且基因之间也存在overlap, 在对这些特殊情况进行处理时,htseq-count内置了以下3种模式

union

intersection-strict

intersection-nonempty

通过–mode参数指定某种模式,默认值为union。这3种模式在判断一条reads是否属于某个feature时,有不同的判别标准,示意图如下

在BAM文件,包含了比对上的reads和没有比对上的reads, 只有比对上的reads 会用来计数,htseq-count默认会根据mapping的质量值对BAM文件进行过滤,默认值为10, 意味着只有mapping quality > 10的reads才会用来计数,当然可以通过-a参数来修改这个阈值。

能够明确reads属于一个featurer时,比如示意图种的第一种情况,

reads完全是gene_A的一个片段,该feature的计数就加1;

能明确reads不属于一个feature时,称之为no_feature, 比如示意图种的第二种情况,

reads的一部分比对上了gene_A, 在intersection_strict模式下,判定该reads不属于gene_A, 就是no_feature。

当不明确一条reads是否属于某个feature时,通常是由于reads在两个gene的overlap区,比如示意图中的第六和第七种情况,这样的reads被标记为ambigous。

当一条reads比对上了两个feature时,会被标记为alignment_not_unique。

在统计属于某个基因的reads数时,需要重点关注对 ambiguous 和 alignment_not_unique 的reads的处理, 通过–nonunique参数来指定,取值有以下两种

none,all

默认值为none时,这两种reads被忽略掉,不进行任何的计数;取值为all时,对应的所有feature的计数都会加1。

除了–mode和–nonunique两个参数外,还需要关注–stranded参数,这个参数指定文库的类型,默认值为yes, 代表文库为链特异性文库,no代表为非链特异性文库。

对于非链特异性文库文库,在判断一条reads是否属于一个基因时,只需要关注比对位置,默认值为yes, 代表文库为链特异性文库,而链特异性文库还需要关注比对的正负链和基因的正负链是否一致,只有一致时,才会计数。

理解了以上3个参数,就能够正确的使用htseq-count了。对于非链特异性的数据,常规用法如下

htseq-count -f bam -r name -s no -a 10 -t exon -i gene_id -m union –nonunique=none -o htseq.count align.sorted.bam hg19.gtf

在运行速度上,featurecounts比htseq-count快很多倍,而且feature-count不仅支持基因/转录本的定量,也支持exon等单个feature的定量。所以更加推荐使用featurecounts来定量。

#htseq-count 命令参数

-f | –format default: sam 设置输入文件的格式,该参数的值可以是sam或bam。

-r | –order default: name 设置sam或bam文件的排序方式,该参数的值可以是name或pos。前者表示按read名进行排序,后者表示按比对的参考基因组位置进行排序。若测序数据是双末端测序,当输入sam/bam文件是按pos方式排序的时候,两端reads的比对结果在sam/bam文件中一般不是紧邻的两行,程序会将reads对的第一个比对结果放入内存,直到读取到另一端read的比对结果。因此,选择pos可能会导致程序使用较多的内存,它也适合于未排序的sam/bam文件。而pos排序则表示程序认为双末端测序的reads比对结果在紧邻的两行上,也适合于单端测序的比对结果。很多其它表达量分析软件要求输入的sam/bam文件是按pos排序的,但HTSeq推荐使用name排序,且一般比对软件的默认输出结果也是按name进行排序的。

-s | –stranded default: yes 设置是否是链特异性测序。该参数的值可以是yes,no或reverse。no表示非链特异性测序;若是单端测序,yes表示read比对到了基因的正义链上;若是双末端测序,yes表示read1比对到了基因正义链上,read2比对到基因负义链上;reverse表示双末端测序情况下与yes值相反的结果。根据说明文件的理解,一般情况下双末端链特异性测序,该参数的值应该选择reverse(本人暂时没有测试该参数)。

-a | –a default: 10 忽略比对质量低于此值的比对结果。在0.5.4版本以前该参数默认值是0。

-t | –type default: exon 程序会对该指定的feature(gtf/gff文件第三列)进行表达量计算,而gtf/gff文件中其它的feature都会被忽略。

-i | –idattr default: gene_id 设置feature ID是由gtf/gff文件第9列那个标签决定的;若gtf/gff文件多行具有相同的feature ID,则它们来自同一个feature,程序会计算这些features的表达量之和赋给相应的feature ID。

-m | –mode default: union 设置表达量计算模式。该参数的值可以有union, intersection-strict and intersection-nonempty。这三种模式的选择请见上面对这3种模式的示意图。从图中可知,对于原核生物,推荐使用intersection-strict模式;对于真核生物,推荐使用union模式。

-o | –samout 输出一个sam文件,该sam文件的比对结果中多了一个XF标签,表示该read比对到了某个feature上。

-q | –quiet 不输出程序运行的状态信息和警告信息。

-h | –help 输出帮助信息。

结果每个样本输出一个count文件,其中包含了基因名和reads计数;另外,如果你看到文件倒数5行,tail htseq.count.txt 会发现还有几行带文字的

no_feature:比对区域与任何基因都没有重叠 。

ambiguous:比对区域与多个基因都发生重叠

 too_low_aQual:比对质量低于设定阈值(默认是10)

 not_aligned:无法比对上

 alignment_not_unique:比对位置不唯一

对两个软件的结果进行合并

##合并htdeq生成的count文件到matrix.count

cd $MATRIX/htseqperl -lne ‘if ($ARGV=~/(.*).count/){print “$1t$_”}’ *.count | grep -v “_” >matrix.count

##合并featureCounts生成的count文件到matrix_2.count

cd $MATRIX/featureCountsfor i in `seq 15 18`;docut -f 1,7 SRR20382${i}.feature.count | grep -v “^#” > SRR20382${i}.matrixsed -i ‘1d’ SRR20382${i}.matrixdoneperl -lne ‘if ($ARGV=~/(.*).matrix/){print “$1t$_”}’ *.matrix >matrix_2.count

统计一下两个软件的计数之和

#统计featureCounts  

perl -alne ‘$sum += $F[2]; END {print $sum}’ matrix.count

#结果是5882943

#统计htseq-count,结果是786338

6人点赞

更多精彩内容,就在简书APP
“小礼物走一走,来简书关注我”

还没有人赞赏,支持一下

总资产2 (约0.14元)共写了6.1W字获得57个赞共69个粉丝

点击数:0

国家卫健委《医疗美容主诊医师备案培训大纲》正式发布

国家卫健委《医疗美容主诊医师备案培训大纲》正式发布

原创

口腔之美



口腔之美


2020-07-15

       国家卫健委组织制定的《医疗美容主诊医师备案培训大纲》于2020年7月6日正式发布。“大纲”是在国家卫健委组织领导下,国家卫健委医政医管局具体指导和政策把关、审核,中国整形美容协会组织完成,是保证中国整形美容健康发展的重要专业培训大纲。国家卫健委办公厅在通知中指出:“为落实《医疗美容服务管理办法》有关要求,指导各地规范开展医疗美容主诊医师培训工作,进一步提高医疗美容主诊医师备案管理水平,切实保障医疗美容服务质量和安全,国家卫健委组织制定了《医疗美容主诊医师备案培训大纲》。供各地参考使用。现将“大纲”中的美容牙科(口腔美容)部分摘录,供口腔同仁学习参考。

       美容牙科部分由刘洪臣等专家起草、修改、审校,历经数稿,经中国整形美容协会口腔整形美容分会及全国3分之2省份口腔专家代表充分讨论提出意见和建议,并在《口腔颌面修复学》杂志2018年第1期刊登征求意见稿,再由中整协组织专家及常务理事会讨论通过,报国家卫健委批准制定。由于《医疗美容服务管理办法》中关于“美容牙科”的名称尚未修改,按规定仍沿用。培训大纲恢复了口腔颌面美容的内容,特别将以往口腔颌面美容及注射填充美容的内容完善。大纲的发布是对包括口腔整形美容的各专业美容主诊医师备案培训的规范,将有许多工作要做。包括大纲的具体实施方案、培训机构的认定、各地卫生行政管理机构的备案要求等。2011年由当时卫生部医政医管司指导,中华医学会与人民卫生出版社组织,刘洪臣教授任总主编的《美容主诊医师培训大纲与方案》及各专业系列教材也要根据颁布的备案培训大纲做相应修改。       

       各位从事口腔美容的口腔同仁注意按照国家卫健委的要求,在认定的培训机构严格按照规定、标准、规范培训合格、在当地卫生行政部门备案后再开展特定的口腔美容整形项目。还要注意学习国家其他关于医疗美容的法律法规,特别是《医疗美容服务管理办法》、《医疗美容分级管理目录》以及医疗美容准人的相关规定,落实国家8部委关于打击非法医疗美容的通知精神,避免发生不必要的医疗过错或纠纷。也请大家提出意见和建议,使大纲的美容牙科部分全面实施。

       《医疗美容主诊医师备案培训大纲》美容牙科专家组                                           中国整形协会口腔整形美容分会  


美容牙科主诊医师备案培训大纲

(共计48学时)

单元

细目

要点

学时

第一章 绪论

第一节 口腔美容医学的定义

一、口腔美容医学的概念及特点

二、口腔美容医学的起源、发展及其与其他相关学科的关系

1

第二节 口腔美容医学的特点

1

第二章 口腔美容医学的基础理论

第三节 牙齿色彩学

一、牙体硬组织结构与功能的生理学和美学意义

二、牙齿色彩的形成

三、陶瓷修复体的光学特征

四、人工牙的比色和选色

五、前牙审美的视觉规律

2

第四节 口腔颌面部美学

一、颌骨组织的结构与功能

二、颌面软组织的结构与功能

三、牙周组织的结构与功能

四、口腔颌面部美学的客观指标

五、软组织美学

六、微笑的解剖学基础

七、数字美在口腔美容医学中的应用

2

第三章 变色牙的治疗

第五节 牙齿染色、变色的原因及分类

一、外源性:烟、茶、咖啡及有色食物染色

二、内源性:四环素牙、氟斑牙及牙髓坏死造成的牙齿变色

2

第六节 牙齿抛光术

0.5

第七节 牙齿喷砂术

0.5

第八节 洁治技术

一、龈上洁治术

二、龈下刮治术

0.5

第九节牙齿漂白技术

▲(牙齿漂白)

一、诊室内、外漂白方法

二、根管漂白技术

三、家庭漂白

1

第十节 变色牙的修复治疗

一、贴面修复

二、冠修复

0.5

第四章 牙体缺损的美容修复

第十一节 概述

一、牙体缺损的病因

二、牙体缺损的修复原则

1

第十二节 树脂充填术

一、适应证和禁忌证

二、充填技术

0.5

第十三节嵌体

▲(嵌体修复)

一、适应证和禁忌证

二、分类

三、嵌体的美学修复设计

0.5

第十四节贴面修复

▲(瓷面修复)

一、适应证和禁忌证

二、分类

三、贴面的美学修复设计

0.5

第十五节桩核修复

一、适应证和禁忌证

二、分类

三、桩核的美学修复设计

0.5

第十六节全冠修复

一、适应证和禁忌证

二、分类

三、全冠的美学修复设计

0.5

第十七节 部分冠修复

一、适应证和禁忌证

二、分类

三、部分冠的美学修复设计

0.5

第五章 牙周美学治疗

第十八节 牙周疾病与口腔美容的关系

一、牙龈炎

二、牙周炎

三、牙龈色素沉积

1

第十九节龈成形术

▲(牙龈成形术)

一、适应证和禁忌证

二、龈成形术的种类和方法

0.5

第二十节 翻瓣术

一、适应证和禁忌证

二、翻瓣术

0.5

第二十一节 牙周诱导再生

▲(牙周引导组织再生术)

一、适应证和禁忌证

二、牙周诱导再生术

0.5

第二十二节 牙龈色素去除术

一、适应证和禁忌证

二、牙龈色素去除术的种类和方法

0.5

第六章 牙列缺损与缺失的美学修复

第二十三节 固定义齿修复

一、适应证和禁忌证

二、固定义齿的组成及分类

三、固定义齿的美学修复设计

1

第二十四节 可摘局部义齿修复

一、适应证和禁忌证

二、可摘局部义齿的分类

三、可摘局部义齿的美学修复设计

1

第二十五节 全口义齿修复

一、无牙颌的解剖学基础

二、全口义齿人工牙的选择和排列

三、全口义齿美学修复设计

1

第二十六节 附着体义齿

一、适应证和禁忌证

二、附着体义齿的分类

三、附着体义齿的美学修复设计

1

第二十七节 种植义齿修复

一、适应证和禁忌证

二、种植义齿的组成

三、种植义齿的美学修复设计

1

第二十八节 赝附体

一、阻塞器的美学修复设计

二、义眼

三、义鼻

四、义耳

1

第二十九节 临时义齿修复

一、临时可摘局部义齿的设计及制作

二、临时固定义齿的设计及制作

三、临时全口义齿的设计及制作

1

第七章 牙饰技术

第三十节 钻石型水晶贴面

一、适应证与禁忌证

二、钻石型水晶贴面的设计及制作

1

第三十一节 其他类型饰齿与文齿

一、其他类型饰齿

二、文齿

三、饰齿文齿的美学评价

1

第八章 错 畸形的矫治

第三十二节 错畸形的分类

一、错畸形的原因

二、错畸形的分类

1

第三十三节 各类错畸形的矫治

一、安氏I类错畸形的矫治

二、安氏II类错畸形的矫治

三、安氏III类错畸形的矫治

1

第三十四节 常用美容矫治技术

1

第九章 口腔颌面部软组织微整形美容

第三十五节 口腔颌面部组织注射美容技术

1

第三十六节 口腔颌面部组织充填美容技术

1

第十章 口腔颌面部整形美容

第三十七节 口腔颌面部整形美容手术基础

一、治疗计划的制定

二、手术入路及切口设计与闭合方法

三、颌面部瘢痕的处理及预防

2

第三十八节 上下颌前突矫正术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

2

第三十九节 下颌角肥大矫正术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

2

第四十节 颏成形术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

2

第四十一节 唇、腭裂修复术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

五、唇腭裂的序列治疗及功能训练

2

第四十二节 颧骨整形美容

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

2

第四十三节 面部小瘢痕整复术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

1

第四十四节 唇、颊软组织整形美容

一、唇、颊软组织常见美容问题

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

2

第四十五节 面部黑痣手术

一、适应证及禁忌证

二、术前准备

三、手术步骤及要点

四、术后处理及并发证的防治

1

信息来源:国家卫生健康委


《医疗美容主诊医师备案培训大纲》美容牙科起草修订专家组

组长:刘洪臣

成员:李鸿波  张志光  沈国芳  周诺  王成龙

秘书:刘乙颖 


中国整形美容协会口腔整形美容分会委员会

主任委员:刘洪臣

候任主任委员:张志光

副主任委员:房兵 田卫东 沈刚 周彦恒 

                       赵吉宏 唐瞻贵 董福生

常务委员:

马净植  王仁飞  王鹏来  王燕一  牛光良 

田卫东 朱洪水 朱娟芳 任贵云 刘青梅 刘峰 刘强 刘静明 米方林 汤春波 许彪 孙迎春 

孙勇 李冰 李新春 李肇元 杨晓江 沈刚 

张月兰 张东升 张旻 张晓东 陈莉莉 陈曦 

邵龙泉 林松杉 周彦恒 周洪 郑东翔 房兵 

赵文峰 赵吉宏 赵华强 赵红宇 赵志河 赵克 

段义峰 段瑞平 贾俊 徐宝华 徐普 高秀秋 

唐瞻贵 黄远亮 葛林虎 董福生 程涛 曾融生 

谢伟丽 谢晓莉 廖红兵 谭颖徽

委    员:

丁加根 丁志江 于书娟 马净植 马俊青 

马菲菲 马攀 王少海  王仁飞 王旭东 

王来平 王俊成 王莉莉 王桃 王涛 王鹏来 

王燕一 王鲲 牛光良 方明 邓邦莲 邓悦 

左金华 田卫东 冯驭驰 冯兴梅 冯剑颖 

曲卫国 朱小峰  朱东望 朱冰生 朱洪水 

朱娟芳 乔光伟 乔栒柏 任贵云 刘长虹 

刘正彤 刘加强 刘冰 刘宇 刘青梅 刘畅 

刘怡 刘洪臣  刘峰 刘雪梅 刘斌 刘强 

刘静明 刘聪 米方林 汤春波 许彪 许跃 

孙迎春 孙贵峰 孙勇 孙晓梅 孙健 花菲 

杜平功 李玉梅 李永斌 李冰 李志革 李丽雅 

李祎 李秋影 李清 李鸿波 李琦 李斌 李新春 

李静 李肇元 李毅萍 杨晓江 杨晓喻 肖金刚 

肖惠军 吴永生 吴晓霞 吴烨 吴楠 邱宜农 

冷卫东 冷斌 闵安杰 汪振华 沈 刚 张月兰 

张文怡 张文禹 张东 张东升 张庆福 张志光 

张君 张旻 张珂  张相皞 张晓东 张健 张娟 

张萍 张琛 张喆焱 张雷 张蕾 陈良娇 陈奕帆 

陈莉莉 陈晓红 陈润 陈雪峰 陈琰 陈曦 

邵龙泉 林云红 林松杉 林婷 欧阳东 罗峰 

周珊 周彦恒 周洪 郑东翔 房兵 赵卫星 

赵文峰 赵冬 赵吉宏 赵华强 赵红宇 赵志河 

赵克 胡骁颖 段义峰 段曰廷 段瑞平 侯敏 

姜涛 姚源 栗洪师 贾俊 铁朝荣 徐宝华 

徐俊峰 徐普 徐璐璐 高秀秋 高涛 郭航 

唐瞻贵 黄文霞 黄兰 黄达鸿 黄远亮 

黄晓峰 崔颖秋 康非吾 章捍东 随丽娜  彭勃 

葛成 葛林虎 董福生 韩建国 程杰 程涛 

傅进友 焦婷 曾融生 温秀杰 谢伟丽 谢奇 

谢晓莉 谢逸瑞  粱国健 裘松波 解建立  

蔡巧玲  廖光天 廖红兵 廖湘凌 廖楚航 

谭颖徽 颜兴 薛启明 戴丽霞 魏利敏

(以上为第二届口腔整形美容分会名单,部分军队委员已按要求申请退出)


编辑:刘乙颖

审核、校对:刘  洋

※《医疗美容主诊医师备案培训大纲》可点击左下角“阅读原文”自行下载

点击数:0

配制100ml完全成骨诱导液

配制100ml完全成骨诱导液
 
 
 
 
 
Stocking 浓度
分装
取液量
终浓度
取液量
     终浓度         
 
地塞米松
3.925mg
无水乙醇
10ml
1mM
1ml*10
1ul
10nM
10ul
100nM
M:176.12
维生素C
50mg
PBS
10ml
5mg/ml
1ml*10
1ml
50ug/ml
 
 
 
Beta-甘油磷酸钠
3.0611g
PBS
9.1ml
1M
1ml/*10
1ml
10mM
 
 
 
 
 
 
 
 
 
 
 
 
 

点击数:0