利用 OSCA 进行 eQTL 分析
Nathan写于20200807。
OSCA(OmicS-data-based Complex trait Analysis)是杨老师2019年发表的用于分析多组学数据复杂性状的软件。
OSCA主要可以做以下几个事情:
- 计算表观标记(或者基因表达)数据个体之间的表观组(转录组)的亲缘关系
- 基于混合线性模型的方法检测DNA甲基化(或者基因表达)标记与复杂性状之间的关联(使用MOA和MOMENT方法)
- 在混合线性模型(类似于BLUP)中,估算所有甲基化(转录本)标记对表型联合的“效应”
- 基于具有相关协变量的线性回归模型的eQTL / mQTL分析
OSCA软件的下载和说明是在一个独立的网站,并没有发表在GitHub,而且这个网站有多个杨老师开发的软件,比如GCTA等,有兴趣的话可以浏览一下https://cnsgenomics.com/software/osca/#Overview
00准备工作
eQTL不是一个新奇的分析了,早在12年就有专门针对eQTL的R发表了——MatrixEQTL。虽然MatrixEQTL仍然是现在很多文章是用的软件,但是由于它是R写的,并且输入文件也比较繁琐和复杂,所以我们这次是使用OSCA去做eQTL。
这里先讲明eQTL的几个概念,cis-eQTL和trans-eQTL。cis-eQTL就是某个基因的 eQTL 定位到该基因所在的基因组区域,表明可能是该基因本身的差别引起的 mRNA 水平变化;后者是指某个基因的 eQTL 定位到其他基因组区域,表明其他基因的差别控制该基因 mRNA 水平的差异。
用OSCA做eQTL的流程很简单,主要麻烦的地方在输入数据的准备上。准备好数据直接跑osca --eqtl就可以了。
01 Input数据准备
OSCA做eQTL需要两个输入数据,一个是表型数据,这个数据的格式是OSCA独特的BOD格式,另一个就是marker了,这个数据格式官网只说了PLINK二进制格式,其他的格式不知道可不可以。
首先来看一下BOD格式,这个数据格式和PLINK比较相似,也是一个文件名后的三个不同后缀的文件为一个整文件。BOD格式由三个不同内容的格式文件组成,首先是oii格式,这个格式类似于PINK的fam文件,其需要五列信息,family ID,individual ID,paternal ID,maternal ID和性别,其中1为男性,2为女性,0则是代表未知,Missing用"NA"代替。
其次是opi格式文件,这个文件记录的是转录组的信息,这个文件是一个特殊的文件,这个数据的目的是将loci和基因位置相结合。这个文件包含了五列数据,分别是染色体,probe ID(官网讲这个ID可以使一个外显子或者转录本的ID),物理位置(这个就很费解,不管是转录本还是外显子他的位置都是区域,这怎么定义呢),然后是基因的ID和基因的方向。看一个实例就明白了。
现在我越来越怀疑,这个位置信息怎么填写,只有一个position的话,怎么知道是cis还是trans呢。发个邮件问一下好啦。
最后一个文件是bod格式的二进制文件。像plink一样,二进制文件我们一般是不可以直接编辑的,osca也同样给出了如何做bod file的命令。
# compile data in binary format from text format
osca --efile myprofile.txt --methylation-beta --make-bod --out myprofile
- --efile reads a DNA methylation (or gene expression) data file in plain text format.
- --methylation-beta indicates methylation beta values in the file.
- --make-bod saves DNA methylation (or gene expression) data in binary format.
- --out saves data ( or results) in a file.
-
--gene-expression indicates gene expression profiles in the file.
这样的话就需要我们给出一个表观数据或者转录组的数据,我们来看一下官网给出的实例数据。
推荐阅读
-
这个数据库收集了 24,000 个 "梦"!它利用人工智能对你的梦境进行分析和评分,找出梦境与现实之间的关联,帮助你解梦。
-
正负偏差变量 即 d2+、d2- 分别表示决策值中超出和未达到目标值的部分。而 di+、di- 均大于 0 刚性约束和目标约束(柔性目标约束有偏差) 在多目标规划中,>=/<= 在刚性约束中保持不变。当需要将约束条件转换为柔性约束条件时,需要将 >=/<= 更改为 =(因为已经有 d2+、d2- 用来表示正负偏差),并附加上 (+dii-di+) 注意这里是 +di、-di+!之所以是 +di,-di+,是因为需要将目标还原为最接近的原始刚性约束条件 优先级因素和权重因素 对多个目标进行优先排序和优先排序 目标规划的目标函数 是所有偏差变量的加权和。值得注意的是,这个加权和都取最小值。而 di+ 和 dii- 并不一定要出现在每个不同的需求层次中。具体分析需要具体问题具体分析 下面是一个例子: 题目中说设备 B 既要求充分利用,又要求尽可能不加班,那么列出的时间计量表达式即为:min z = P3 (d3- + d3 +) 使用 + 而不是 -d3 + 的原因是:正负偏差不可能同时存在,必须有 di+di=0 (因为判定值不可能同时大于目标值和小于目标值),而前面是 min,所以只要取 + 并让 di+ 和 dii- 都为正值即可。因此,得出以下规则: 最后,给出示例和相应的解法: 问题:某企业生产 A 和 B 两种产品,需要使用 A、B、C 三种设备。下表显示了与工时和设备使用限制有关的产品利润率。问该企业应如何组织生产以实现下列目标? (1) 力争利润目标不低于 1 500 美元; (2) 考虑到市场需求,A、B 两种产品的生产比例应尽量保持在 1:2; (3)设备 A 是贵重设备,严禁超时使用; (4)设备 C 可以适当加班,但要控制;设备 B 要求充分利用,但尽量不加班。 从重要性来看,设备 B 的重要性是设备 C 的三倍。 建立相应的目标规划模型并求解。 解:设企业生产 A、B 两种产品的件数分别为 x1、x2,并建立相应的目标计划模型: 以下为顺序求解法,利用 LINGO 求解: 1 级目标: 模型。 设置。 variable/1..2/:x;! s_con_num/1...4/:g,dplus,dminus;!所需软约束数量(g=dplus=dminus 数量)及相关参数; s_con(s_con_num);! s_con(s_con_num,variable):c;!软约束系数; 结束集 数据。 g=1500 0 16 15. c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(1);!第一个目标函数;!对应于 min=z 的第一小部分;! 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); !使用设置完成的数据构建软约束表达式; ! !软约束表达式 @for(variable:@gin(x)); !将变量约束为整数; ! 结束 此时,第一级目标的最优值为 0,第一级偏差为 0: 第二级目标: !求 dminus(1)=0,然后求解第二级目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=dminus(2)+dplus(2);!第二个目标函数 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标结果 @for(variable:@gin(x)); ! 结束 此时,第二个目标的最优值为 0,偏差为 0: 第三目标 !求 dminus(2)=0,然后求解第三个目标。 模型。 设置。 变量/1..2/:x;!设置:变量/1..2/:x; ! s_con_num/1...4/:g,dplus,dminus;!软约束数量及相关参数; s_con(s_con_num(s_con_num));! s_con(s_con_num,variable):c;! 软约束系数; s_con(s_con_num,variable):c;! 结束集 数据。 g=1500 0 16 15; c=200 300 2 -1 4 0 0 5; 结束数据 min=3*dminus(3)+3*dplus(3)+dminus(4);!第三个目标函数。 2*x(1)+2*x(2)<12;!硬约束 @for(s_con_num(i):@sum(variable(j):c(i,j)*x(j))+dminus(i)-dplus(i)=g(i)); ! 软约束表达式;! dminus(1)=0; !第一个目标约束条件; ! dminus(2)+dplus(2)=0; !第二个目标约束条件 @for(variable:@gin(x));! 结束 最终结果为 x1=2,x2=4,dplus(1)=100,最优利润为
-
利用 OSCA 进行 eQTL 分析
-
[Matlab 量化投资] 利用数据包络分析和遗传算法进行选股分析?你懂了吗?(附源程序)
-
如何利用数据包络分析法(DEA)进行效率评估?
-
仪表盘显示 | 利用 DataEase 对 2022 年北京冬奥会数据进行可视化分析
-
利用 ArcGIS + Fragstats 软件进行景观模式指数分析(附练习数据下载)
-
GIS 云平台 利用云计算进行 GIS 效益分析
-
第 20 天:NLP 实践(四)--利用 GRU 模型对电影评论进行情感分析
-
利用 SWOT 分析和蓝海战略进行营销巡演