生信专用简明 Python 文字和视频教程 源码在:https://github.com/Tong-Chen/Bioinfo_course_python Reference 一些练习题 给定FASTA
导言 构建生信分析流程是生物信息学从业人员必备的技能之一,对该项能力的评估常常是各大公司招录人员的参考项目之一。 一个好的生物信息分析流程可以让你事倍功半,有效减负,同时也有利于他人重复你的数据分析结果。 ? 根据生信信息学数据分析流程(管道、工作流程序)构建的风格和方式,大致有以下几大流派(注1): 脚本语言流 Common Workflow language 语言流 Makefile流 配置文件流 Jupyter 生信分析流程构建的几大流派 | 脚本语言流 脚本语言流的主要是通过简单的脚本语言(如shell,R,Python,Perl)运行各类命令行脚本/程序。 )最早开始接触生信分析流程的方式。
导言 构建生信分析流程是生物信息学从业人员必备的技能之一,对该项能力的评估常常是各大公司招录人员的参考项目之一。 一个好的生物信息分析流程可以让你事倍功半,有效减负,同时也有利于他人重复你的数据分析结果。 根据生信信息学数据分析流程(管道、工作流程序)构建的风格和方式,大致有以下几大流派(注1): 脚本语言流 Common Workflow language 语言流 Makefile流 配置文件流 Jupyter 生信分析流程构建的几大流派 | 脚本语言流 脚本语言流的主要是通过简单的脚本语言(如 shell,R,Python,Perl)运行各类命令行脚本/程序。 前两种(1 和 2)是大多数生物信息学初学者(不具备封装和打包能力)最早开始接触生信分析流程的方式。后两种(3 和 4)是专业人员开发新工具、新流程的必备技能。
接上两篇内容,本文主要讲述工作中NGS从科研进入医学临床领域,工作中接触到生信流程,以及最终在实现的过程。 接触二代测序,生信分析,那真是打开了一个新世界的大门,各种名次术语满天飞,搞的头晕脑胀。 下面分阶段描述生信分析流程升级/进化的过程: ---- 1.手动命令行运行 经过几个月接触,自学、爬坑,慢慢搞清楚了部分内容,在似懂非懂之间开始了生信流程分析,终于有一天明白过来,这所谓的pipeline 比如其中一个步骤: 生信分析流程的进化_1.png QC 完成后,然后运行下一个步骤: 生信分析流程的进化_2.png 运行模式,一个输入或者多个输入文件,通过软件分析/计算得到一个或者多个输出文件 脚本连续运行 随着熟练程度提高,生信分析上用到的软件/工具也熟悉起来了,但是问题也暴露出来了,简单的一套 GATK Best Practice 肿瘤突变分析流程,加上CNV,SV 分析从 fastq 文件开始到最后得到过滤的 到这里,基本上就达到很多公司的生信自动化分析水平了 6. 然而到这里就足够了么?
通过分析 CDR3 motif 的序列特征(如长度、关键氨基酸残基),可阐明免疫受体(Ig/TCR)与抗原结合的分子规律,揭示 “抗原 - 受体” 相互作用的特异性机制,为理解免疫应答的精准调控提供依据 图1.不同长度的CDR3 氨基酸motif比较 2.2 CDR3 中间氨基酸motif 分析 这是一种忽略CDR3长度,关注CDR3中间氨基酸多样性motif分析策略。 图2.CDR3中间氨基酸motif 分析 3. 举例分析CDR3motif 3.1 展示特定长度CDR3 氨基酸motif library(tidyr) library(tidyverse) files <-list.files(path = ". page_grobs) + 1]] <- combined_plot } # 将所有页面组合成一个 PDF 文件 pdf(paste0(condition_name,".png"), width = 9,
第一步,下载COAD数据 ########################################################################################## ## step1 load package and change Working Directory ########################################################################################### library(
在上文图形化开放式生信分析系统开发 - 3 生信分析流程的进化 讨论了生信分析pipeline的进化,从手动到自动,但仍然停留在终端命令行阶段,为了让更多非生信专业的人能够使用,就要想办法实现生信分析pipeline 变量的分类: 根据实际经验,pipeline变量用到根据用途可以分为以下几类: 分析数据目录${data} 分析过程输出目录${result} 分析用到的软件 bwa ${tools.bwa};samtools ${tools.samtools} 分析流程中用的reference文件以及数据库,如 hg19.fa ${ref.hg19} 分析流程中,用到的cutoff值. 如 cnv的cutoff值 ${cutoff.cnv} 分析流程运行时配置的资源,如 线程数 ${threads} 分配内存大小${mem} 变量值的类型: 字符:通用的格式,比较宽松 /节点设计:前文提到,生信分析pipeline其实就是基于文件输入输出的工作流,这里对工作流做了简化,归纳起来工作流中有4种节点。
本篇内容引自生信技能树 DAY7-9 课前提问: 1、为什么要做数据挖掘? 即用别人的数据用在自己的文章里面,多半是从别人的数据里筛选自己想要的基因。 拿到表达矩阵之前的分析叫上游分析,之后的叫下游分析。 5、筛选基因的分析方法(基因芯片和转录组通用) 差异分析:测定基因表达量的原理不同,前者是光信号值,后者是短序列的条数。 三、代码分析流程(以二分组代码举例) 1、下载数据,提取表达矩阵和临床信息 rm(list = ls()) #打破下载时间的限制,改前60秒,改后10w秒 options(timeout = 100000 只是芯片能用,其它组学用不了 exp <- exprs(eSet) #⭐第一个要检查的地方,表达矩阵行列数,正常是几万行,列数=样本数, #行代表探针数,如果0行说明不是表达芯片或者是遇到特殊情况,不能用此流程分析 step2output.Rdata") #比较复杂的探针注释参考资料 #资料1:拆分取列https://www.yuque.com/xiaojiewanglezenmofenshen/kzgwzl/sv262capcgg9o8s5
差异分析表格二分组数据差异分析#差异分析 limmalibrary(limma)design = model.matrix(~Group) # 生成模型矩阵fit = lmFit(exp,design) ,]rownames(exp3) = ids$symbolexp3[1:4,1:4]exp4 = limma::avereps(exp3)# 此时拿到的exp4已经是一个基因为行名的表达矩阵,直接差异分析 4.加ENTREZID列,用于富集分析(symbol转entrezid,然后inner_join)library(clusterProfiler)library(org.Hs.eg.db)s2e = bitr nrow(deg)save(exp,Group,deg,logFC_t,p_t,file = "step4output.Rdata")富集分析# 准备工作rm(list = ls()) load(file clusterProfiler-book/index.html# GOplot:https://mp.weixin.qq.com/s/LonwdDhDn8iFUfxqSJ2Wew# 网上的资料和宝藏无穷无尽,学好R语言慢慢发掘~生信技能树
生信论文的套路 ONCOMINE从全景、亚型两个维度做表达差异分析; 临床标本从蛋白水平确认(或HPA数据库),很重要; Kaplan-Meier Plotter从临床意义的角度阐明其重要性; cBio-portal 数据库做基因组学的分析(机制一); STRING互作和GO/KEGG分析探讨可能的信号通路(机制二); TISIDB/TIMER分析肿瘤免疫特征(机制三)。 差异分析,无论是Oncomine,GEPIA,还是UALCAN、HPA数据库,都不需要R语言编写代码,容易上手,基本上一个星期甚至更短的时间就可以搞定,属于菜鸟级别生信操作。并没有想象中那么难。 生存分析是生信论文中经常出现的表型,也就是说基因在正常和肿瘤组织中表达的差异,与生存率的指标密切相关。如A基因在肿瘤中表达明显上调,生存率显著下降,这就是非常明确的相关性。 生存分析是非常重要的表型,诸多文章均有介绍。这里,我们对生存分析的纯生信数据库进行总结,果友们在选择时也可以作为参考。
好在时至今日,已经有很多科学家开发了非常多优秀的算法及软件,很多时候我们要做的是怎么将这些软件串联起来并构建成生信分析流程,而这项技能通常是各大公司考核应聘者的项目之一。 这种理念是我们最直观的分析逻辑,也是最常用的流程框架。通常,生信刚入门的同学们会选择这种方式,简单而暴力;段位较高的同学,则会选择将分析内容进行包装,然后提供多个参数选择,增加流程灵活性。 大部分时候,这样都会满足我们分析需求,但是其作为一个生信流程有着严重的缺点就是缺乏重入性(reentrancy),即当流程在运行过程中,很容易因为某些不知名的原因而发生中断,而普通的脚本流程只能是从头来过了 下面是Galaxy在线编辑WES分析流程界面: ? ,那么就可以使用Implicit/Explicit类的流程,如:Snakemake、Nextflow等,而这一类的流程也比较适合刚入门生信的小伙伴们去尝试; 如果是需要进行高性能流程开发,致力于解决特定的生物学问题
在上文图形化开放式生信分析系统开发 - 4 生信分析流程的图形化设计 讨论了生信分析pipeline的图形化,如何用图形的方式显示生信pipeline,但是pipeline脚本按照变量的形式保存之后,如何运行 本程序作为控制端,可以和分析端部署在一台机器,也可以通过联网方式连接。 分析过程状态、服务器运行状态,需要由服务器端推送到用户端。 SliverWorkspace_V2.1.pdf 首先这里实现了,服务器账户信息的管理,账户、主机名、端口、密钥、密码,这些信息为了保证安全,需要二次加密,不能将密码明文保存在数据库中,一旦泄漏危害巨大 服务器信息2.png 针对分析流程 运行完成后服务器端推送信息到控制端,判断是否符合要求,输出文件是否存在 运行失败后服务器端推送信息到控制端,显示错误信息,错误日志,便于生信开发人员查找错误 统计每一个分析步骤的运行时间,便于统计分析
在上文生信分析系统开发 - 5 生信分析流程服务器端运行 解决了设计好的流程在分析服务器上运行的问题,随之而来就衍生出的新需求: 一、批量分析流程的运行顺序 流程输入文件是按照样本编号来匹配运行的,运行顺序就取决于样本信息 样本优先级:此处设计样本录入时候,输入优先级字段范围1-9,数字越大优先级越高 样本顺序:样本默认按照录入的时间先后排序,相同优先级的样本数据,按照录入先后顺序排序运行 二、分析流程的过程控制 强制停止 :开始分析流程之后,可以在任意时间停止分析过程 错误恢复运行:分析流程运行错误,再次运行,可以选择从错误处恢复运行 调试恢复运行:分析结束,修改pipeline参数后,可以选择从任意一点开始重新分析(应用修改后的参数 ) 覆盖运行:分析错误或者正常结束后,重新从头开始运行整个分析流程,并覆盖之前的分析结果 技术实现: 一、批量样本分析流程的运行顺序: 依赖于样本的优先级和输入的顺序先后,每次从数据库取出未完成分析的样本列表 每次列表中取最上面一个,去指定目录里匹配输入文件,如果输入文件符合匹配要求,立即启动pipeline,开始流程分析。同时更新该样本状态为正在分析,这样下次排序时候就不会重复获取该样本了。
参数data:一个数据框,用于做生存分析的数据。 参数p.adjust.method:p值矫正方法(参见p.adjust)。 0.062 0.485 0.527 P value adjustment method: BH 通过这个矩阵我们可以看到molecular_group 水平上两两之间的生存分析对比的
生信论文的套路 ONCOMINE从全景、亚型两个维度做表达差异分析; 临床标本从蛋白水平确认(或HPA数据库),很重要; Kaplan-Meier Plotter从临床意义的角度阐明其重要性; cBio-portal 在差异分析的前提下,表型分析成为重点内容,也是可以玩出花样的地方。 生存分析是非常常见的表型分析。与生存分析相比,相关性分析是另外一个常见的表型分析。 免疫浸润分析比生存分析、差异分析和相关性分析难度更大,因为免疫学是不断延伸、拓展的学科,并不断从理论走向应用、临床,兼有科学性和技术性,比如实验中已经普遍应用的免疫印迹(WB),流式分析,免疫组化和免疫荧光等等 生信分析中,有一种算法叫反卷积分析,英文名叫Deconvolution。 生信开发人员可以先通过预设一个优秀的数据训练集(训练集主要包含了每种不同免疫细胞的基因表达特征),然后通过反卷积算法推算出这个整体样本中究竟有哪些免疫细胞。
上课的好处就是可以听到很多细节和经验,是自己零散地找资料不能相比的,收获很多,感觉要全部吞下来还要再复习几遍 贴具体的代码好像挺重复的,只贴一些过去自学没有了解到的有趣的TIPS: 1、转录组数据不能直接套芯片的流程 ,错套的第一个表现就是dim结果是 行 为0 2、non-coding和普通array可以统一处理,但不能做富集分析,富集分析需要用编码蛋白做,或者先靶基因预测然后再做富集分析 3、normalize不是一定要加 里的第二栏里,带“--”说明不对应任何symbol,需要删去 7、一个探针对应多个基因(非特异性探针),难以解释,这些行直接去掉 8、对于lnkRNA不能直接用页面上的TargetID,尽量寻找symbol列 9、 ="";table(k2) ids = b[k1&k2,] 16、筛选下载数据中的部分样本进行数据分析 library(stringr) # 方法1:按照行号,能数的时候可以自己数行号 keep = c str_detect(pd$title,"and");table(keep) exp = exp[,keep] pd = pd[keep,] 总结: 小洁老师讲的真的非常好,细腻又形象,准备学生信的一定要来学一波让你丝滑入门
每月一生信流程栏目灵感来自于《铁汉1991》博客的《每日一生信》,他那个时候介绍的主要是生信基础知识,包括数据结构,数据格式,数据库资源,计算机基础等等,所以每天都可以进步,每天都有成果。 这些基础知识已经被分享的七七八八了,所以我这里推陈出新,来一个每月一生信流程,陪生信技能树的粉丝们一起进步! image.png 不同数据变换公式的差异 学习这样的流程是需要一定背景知识的 首先是LINUX学习 我在《生信分析人员如何系统入门Linux(2019更新版)》把Linux的学习过程分成6个阶段 , 30万学习量的基础合辑: image.png 生信技能树关于RNA-seq上下游数据分析的教程的确不少了 因为做目录确实很浪费时间,差不多就下面这些,大家先学习吧: 转录组经典表达量矩阵下游分析大全 可变剪切 生信小白的RNA-seq实战历程 RNA-seq数据分析指南 后记 听说隔壁openbiox团队在组织翻译这个bioconductor流程系列,而且还是由我们生信技能树元老-思考问题的熊领头,希望他们的翻译成果早日出版
单基因生信分析流程(1)一文解决TCGA数据下载整理问题 单基因生信分析流程(2)一文解决差异分析和基因相关分析问题 本文目的 学会如何使用差异分析 学会绘制火山图和热图 学会如何求取相关基因 第一招: 差异分析 差异分析步骤总结 (1)读取基因表达矩阵 (2)根据基因表达量设置样本分组 (3)设置差异倍数、生成差异分析结果 (4)绘制火山图和热图 加载所必须的包 # ============== ======================= 设置分组,我们根据ERBB2基因的表达中位值,将样本分为ERBB2高表达组和ERBB2低表达组,通过求两组样本的差异基因,来对ERBB2的生物学功能进行分析
(唯一不足是Conda安装的singularity不是最新版本) conda create -n singularity singularity=3.8.6 -y 4生信分析中为什么使用Singularity Singularity 容器可以将软件及其依赖项打包在一起,确保分析环境的一致性,无论是在本地机器、服务器还是云平台上。 这种高度的移植性确保了生物信息学分析可以轻松地从一个计算环境迁移到另一个,无论是从个人电脑到云平台,还是在不同的研究机构之间共享。 is a test file." > testfile ##创建定义文件 vim test.def ##定义文件写入以下内容 bootstrap: docker from: almalinux:9 singularity build --fakeroot test.sif test.def ## 运行 singularity run test.sif 示例演示 运行输出 6来看一个实例 刚好前段时间生信技能树有发一个推文介绍了中国人群的肝癌多组学队列研究
生信分析,无论是Oncomine,GEPIA,还是KM Plotter数据库,都不需要R语言编写代码,容易上手,基本上一个星期甚至更短的时间就可以搞定,属于菜鸟级别生信操作。并没有想象中那么难。 p<0.05就是我们前期做分析想要的结果,即使用GEO数据库、TCGA数据库做分析,甚至做芯片或测序,没有差异也是枉然。 对于医生,个人建议是最好学会R语言,最起码掌握不用R语言的数据库分析。 最重要的是阅读过的文献和做过的生信分析,可以进一步促进对临床上疾病的认识,提升医生的诊治水平。 ? 接下来,我们将按照中心法则和生信论文分析的思路总结生信分析的网站。 这些网站均是纯生信数据库,不需要R语言基础,就像Excel、PPT一样容易上手,只要你愿意学,肯定能学会。 差异分析数据库 oncomine数据库(差异分析首选) https://www.oncomine.org/resource/main.html GEPIA数据库(共表达是特色) http://gepia.cancer-pku.cn