利用PlinK、EMMAX、R做全基因组关联分析与绘制曼哈顿图、qq图全流程

一、理论部分

二、实操

表型处理

使用Phenotype包(R package)来进行BLUP计算得到校正后的表型值

https://cloud.tencent.com/developer/article/1677209

EMMAX软件

EMMAX的表型数据格式

第一列,FAMID通常指的是“Family ID”的缩写
第二列,INDID代表的是个体的标识符,它是“Individual ID”的缩写,也就是他允许有重复
第三列,表型值

EMMAX的基因组数据格式

vcf文件
前几行是以"#"开头的注释信息,每行是一个SNP,每列是一个个体。行与列整合了一个群体的SNP全部变异信息。
利用plink把vcf文件转换为tped文件。

TPED文件和PED文件的主要区别在于数据的组织和展示方式。以下是具体的比较:

数据组织方式:TPED文件是把PED文件中的信息重新整合,具体做法是在MAP文件的四列之后,加上了PED文件中后面的基因型信息旋转90°之后的形式。因此,TPED文件的一行代表一个SNP的信息,而PED文件的一行代表一个样本的信息。
数据展示方式:在TPED文件中,每行代表一个SNP位点,列数取决于样本量大小,前四列内容同MAP文件,后面所有列为所有样本在该SNPs位点处的基因型信息。而在PED文件中,每行代表一个样本,前六列的内容与TFAM文件相同,后面的列为该样本在各个SNP位点处的基因型信息。
ped文件包含以下几列:

第一列:Family ID。

第二列:Individual ID。自然群体这列和Family ID是一样的。

第三列:Paternal ID。未提供信息的话这列为0。

第四列:Maternal ID。未提供信息的话这列为0。

第五列:Sex。未提供信息的话这列为0。

第六列:Phenotype。一般来说,直接拿vcf转换的话这列为-9,也就是缺失。

第七列开始就是个体在每个标记位点的基因型。

在这里插入图片描述

map文件包含以下几列:

第一列:染色体编号。

第二列:SNP编号。

第三列:遗传距离。未提供信息的话这列为0。

第四列:物理位置。

数据格式与其相互转换

vcf格式转换为ped格式
./plink --vcf /mnt/g/bioinfor/process/GWAS/基因型相关文件/517_plink.vcf --recode 12 --allow-extra-chr --chr-set 26 --out /mnt/g/bioinfor/process/GWAS/基因型相关文件/517

./plink是由于plink软件就在我所在的目录内,若在其他目录,需要指定路径。
–recode参数,设定输出的为为文本格式ATCG,而非二进制。–recode后还可以跟参数,设置SNP单倍型的表现方式,如双数字0/1、0/0或单数字0,1,2,具体见此链接的测试:https://www.cnblogs.com/liujiaxin2018/p/13795230.html
命令运行结束后会输出.map与.ped文件
我用了–recode 12 参数
–chr-set 26设置染色体数目
–allow-extra-chr 允许scaffold染色体

https://blog.youkuaiyun.com/Taylent/article/details/102523295

从ped格式文件中筛选出所需的材料
./plink --file /mnt/g/bioinfor/process/GWAS/基因型相关文件/517N --keep /mnt/g/bioinfor/process/GWAS/基因型相关文件/select.txt  --recode 12 --out /mnt/g/bioinfor/process/GWAS/基因型相关文件/377

–keep功能 需要准备原始的被提取文件,然后是提取需要的参考文件,储存了需要提取哪些样本的信息
这里面用作筛选的select文件的保存格式需要注意,从表格导出txt文本的时候有以制表符分割和无格式文本两种,制表符分割的文件要大一些,是可以工作的,无格式文本会报错:无法读取select文件。
最后是使用PRN格式的文件成功的。
另外第一列是材料的编号,第二列也是材料的编号,只有一行是不行的。这两列在人类或哺乳动物的研究中是用来分家系的,两列可以不同,在植物学的自然群体中,设为相同。
参考:https://dengfei2013.gitee.io/plink-cookbook/plink%E6%96%87%E4%BB%B6%E6%8F%90%E5%8F%96.html#%E6%A0%B7%E6%9C%AC%E6%8F%90%E5%8F%96

ped格式转化为tped格式
 ./plink --file /mnt/g/bioinfor/process/GWAS/基因型相关文件/375 --recode 12 transpose --out t377 --allow-extra-chr --chr-set 26

–file 参数 只需要用户指定一个ped map 文件共有的文件前缀,它会自动读取两个文件。如果按tab补全文件名,则会出错。
–file参数可以替换为 --ped 只转换一个ped文件,这个时候补全ped文件的全名。
不需要指定输出文件的路径,默认为ped文件所在的位置

构建亲缘矩阵
./emmax-kin-intel64  /mnt/g/bioinfor/process/GWAS/基因型相关文件/t375 -v -d 10 -o t377-kin

这里输出文件的位置需要注意,我重复计算了很多遍,因为我找不到输出文件,软件说的-o参数可以指定输出路径与文件名,我发现不管怎么设置,输出文件最后都存放在./emmax-kin-intel64所处的目录,还好计算成功了。

再进一步调整数据

把不同文件的每一行都对应起来。

关联分析计算

/emmax-intel64 -v -d 10 \ #
-t snp_filter \ # 输出snp数据的前缀
-p sample.table \ # 表型数据
-k pop.kinship \ # 亲缘关系矩阵
-o emmax.out \ # 输出文件前缀

./emmax-intel64 -v -d 10 -t /mnt/g/bioinfor/process/GWAS/基因型相关文件/t375 -p /mnt/g/bioinfor/process/GWAS/表型文件/375_adjmean.txt -k /mnt/g/bioinfor/process/GWAS/基因型相关文件/t375-kin -o EMMAX_out

EMMAX 输出计算结果

整理输出结果准备绘图

EMMAX的输出文件后缀为.ps,文件是没有列名的,但是已知第一列为ID,第二列为回归系数(beta),第三列为回归系数的标准差,第四列为为P值。
由于采用的是CMplot包绘图,这个绘图需要构建一个新文件,它的格式是第一列为ID,第二列为chr,第三列position,第四列为p值。第一列的ID与染色体编号都可以从最初从517份重测序文件筛选出的375份材料的map文件,map文件第一列是染色体编号,第二列是SNP的ID,第四列是SNP的物理位置。

map_375<-read.table("G:/bioinfor/process/GWAS/基因型相关文件/375.map")#读入R
plot<-cbind(map_375[,2],map_375[,1],map_375[,4],EMMAX[,4])#重新组合出新数据
plotGWAS<-as.data.frame(plot)#把matrix数据格式转换为data frame

曼哈顿图的阈值计算
在这里插入图片描述
CMplot包自动计算曼哈顿图上显示的阈值,只需提供1/n即可。

CMplot(plotGWAS, threshold = threshold, amplify = F, file = "tiff", plot.type=c("m","q"))

分别绘制曼哈顿图和QQ图,存放到R的工作空间。

### 关于PAT乙级1031题的Python解法 对于PAT乙级中的第1031题,虽然未直接提及该具体题目编号的内容,可以借鉴其他相似难度和类型的题目解答思路来推测解决方案。在处理涉及数据结构操作或者算法设计的问题时,选择合适的数据结构至关重要。 针对链表反转这一类问题,在C++中可以通过指针灵活操控节点完成任务[^1];而在Python里,则可利用列表(list)模拟链表的行为并简化语法表达。尽管存在性能上的差异——由于Python作为解释型语言在执行含有大量循环或I/O操作的任务时效率较低,但这并不妨碍通过优化算法逻辑达成目的。 考虑到实际应用场景下的需求,如果追求更高效的实现方式,也可以考虑采用`collections.deque`这样的双端队列容器类型来进行特定功能开发[^3]。下面提供一段基于Python编写的简单示例代码用于演示如何高效地解决类似问题: ```python from collections import deque def reverse_list(head): node = head stack = [] while node is not None: stack.append(node.value) node = node.next reversed_head = ListNode(-1) current_node = reversed_head while stack: value = stack.pop() new_node = ListNode(value) current_node.next = new_node current_node = new_node return reversed_head.next class ListNode: def __init__(self, val=0, next=None): self.val = val self.next = next ``` 此段代码展示了怎样借助栈(stack)辅助结构逆序遍历单向链表,并重新构建新的反向链接关系。需要注意的是这里为了便于理解采用了较为直观的方式定义了`ListNode`类表示链表结点。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值