【教程】在Radmixture运行无标杆成分名的任意祖源计算器的方法

基因检测 同时被 3 个专栏收录
5 篇文章 0 订阅
2 篇文章 0 订阅
1 篇文章 0 订阅

本文介绍一个使用radmixture程序计算个人祖源方法,如果您做过芯片级基因检测并下载了原始数据(raw data)文件,这篇教程将会给你带来很大的帮助。


一、准备工作

1. 下载自己的基因原始数据。如果您做的是普通的芯片测序,这需要从对应测序机构的网站提供的下载路径与方法的查找;如果您做的是全基因组测序,则先提取BAM/CRAM文件的数据到TXT/CSV的芯片格式里,再进行后续操作即可。

下载、安装R语言。建议从官网下载安装包:https://www.r-project.org/

2. 下载并解压Admixture Studio程序:

https://dnagenics.com/admixture-studio/

部分电脑可能不支持该软件,这需要安装缺少的系统文件来使程序顺利打开。本教程以Lukasz Macuga的LM Genetics K47祖源计算器为例(后文将直接简称为K47)。打开Admixture Studio后在左下方的列表框中选择一个计算器,点击RUN下载并运行计算器。这会同时运行Admixture Studio里的祖源计算器并产生结果,但是这不属于此教程的范围,所以您也可以在计算器下载完后就直接关闭Admixture Studio,主要内容先看后文。(注:如果您能够从其它来源下载到祖源计算器文件的话,也可以直接用而不通过Admixture Studio下载)

3. 接下来在Office Excel或WPS表格中操作(如果您掌握了编程,也可以用更好的方法来替代此方法):

(1) 打开计算器文件所在文件夹,并复制一份计算器的“.alleles”文件;

(2) 修改文件后缀名为Excel或WPS表格格式,打开后按空格分列,然后分别在D1和E1格处输入如下的Excel函数IF,分别在第4、5列互补出第2、3列的基因型:

=IF(B1="A","T",IF(B1="C","G",IF(B1="G","C",IF(B1="T","A",0))))
=IF(C1="A","T",IF(C1="C","G",IF(C1="G","C",IF(C1="T","A",0))))

(3) 最后向下填充、保存。 


二、正式运行radmixture

接下来的步骤需要用到R语言,几乎没有编程基础的用户也可以学习以下内容,引用时需要根据自己的文件位置、祖源计算器信息来替换语句中的彩色加粗字体部分。

1. 安装radmixture,建议从网络端下载;(初次使用只需安装一次,再次调用只需从第2步开始)

install.packages("radmixture")

2. 载入radmixture;

library(radmixture)

3. 读取预备工作中改好的计算器位点文件“.alleles”,注意文件最后写的是什么缀名;

alleles<-read.table(file='X:\\某文件路径\\test.47 - 副本.alleles.xls')

4. 读取计算器文件的频率“.F”(或者“.P”),并输出到中间变量(比如temp1),然后由读取的计算器频率文件生成计算矩阵,接着移除用完的变量;

temp1<-read.table(file='X:\\某文件路径\\test.47.F')
freq<-as.matrix(temp1)
remove(temp1)

5. 读入芯片格式的原始数据,这里需要根据文件格式分两种情况,选择其中一种即可;

(1) TXT格式

genotype<-read.table(file='X:/某路径/某原始数据.txt')

(2) CSV格式

genotype<-read.csv(file='X:/某路径/某原始数据.csv')

6. 将芯片数据对齐到计算器频率矩阵,需要写对祖源成分种类数(以K47为例);

res<-tfrdpub(genotype,47,alleles,freq)

7. 开始频谱分析(由于输出的结果无标杆成分名,pubdata也可以省略);

ances<-fFixQN(res$g,res$q,res$f,tol=1e-5,method="BR",pubdata="")

8. 显示不含标杆成分名的祖源计算器数值结果;

ances$q

如果想直接输出百分比值(不含百分号),则需输入:

100*ances$q

 9. 得到结果后,把成分数值复制到表格文件里,并从对应的计算器文件复制出标杆名,以对应radmixture跑出来的祖源成分数值。然后设计成统计图,设计样式可根据审美要求来决定。


三、拓展内容 

1. 用R语言读取文件路径时,分隔符用“/”或“\\”符号皆可;

2. 如果需要像DIY Dodecad 2.1程序那样读取数据的基因位点使用情况,则需要增加以下的程序语句:

(1) 找出计算器对原始数据的位点利用总数

ncol(ances$f)

(2) 读取计算器的位点总数

nrow(alleles)

(3) 计算位点利用率(genotype rate)

ncol(ances$f)/nrow(alleles)

或者输出为不含百分号的百分数形式:

100*ncol(ances$f)/nrow(alleles)

3. 若要从本地安装已经下载好的radmixture包,则需输入:

install.packages("X:\\某路径\\radmixture_0.0.1.zip")

4. 如果需要卸载ramixture包,则需输入:

remove.packages("radmixture")

5. 如果需要在同一个工作目录文件夹下使用for语句批量计算数据的计算器祖源成分(此时工作目录文件夹内不要放芯片数据以外的文件),则可以参考并感悟以下使用Eurogenes K36祖源计算器做批量统计的示例语句:

library(radmixture)

alleles<-read.table(file='E:\\Admixture\\AdmixtureStudioV250\\Calculators\\EUROGENES_K36_V1\\admix_NEW.alleles')

Temp<-read.table(file='E:\\Admixture\\AdmixtureStudioV250\\Calculators\\EUROGENES_K36_V1\\admix.36.F')
freq<-as.matrix(Temp)
remove(Temp)

setwd('F:/遗传类/现代数据/新建文件夹')
fileNames <-dir()
for (v in fileNames[]){
genotype <- read.table(file = v)
res <- tfrdpub(genotype,36,alleles,freq)
ances <- fFixQN(res$g, res$q, res$f, tol = 1e-5, method = "BR")
results <- cbind(v,ncol(ances$f)/nrow(alleles),100*ances$q)
write.table (results, file ="C:\\admin\\desktop\\K36批量统计结果.txt", sep ="\t", row.names =F, col.names =F, quote =F, append=T)
}

其中setwd语句所在行为设置工作目录,for循环内的write语句所在行为输出计算结果的文本文档目录。


本文完。

本文首发于知乎,同人作品也将同步发表在微基因社区、源基因论坛、CSDN及其他平台。

特别感谢radmixture程序包的出处:https://github.com/wegene-llc/radmixture

  • 0
    点赞
  • 0
    评论
  • 0
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 游动-白 设计师:白松林 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值