标签: genetics

我有一个300 GB的文本文件,其中包含超过250k记录的基因组数据.有些记录包含不良数据,我们的基因组程序'Popoolution'允许我们用星号注释掉"坏"记录.我们的问题是我们找不到将加载数据的文本编辑器,以便我们可以注释掉不良记录.有什么建议？我们有Windows和Linux盒子.

更新:更多信息

Popoolution程序(https://code.google.com/p/popoolation/)在达到"错误"记录时崩溃,向我们提供我们可以注释掉的行号.具体来说,我们从Perl收到一条消息,上面写着"F#€%&Scaffolding".手册建议我们可以使用星号来注释坏线.可悲的是,我们必须多次重复这个过程......

还有一个想法......是否有一种方法可以让我们在不打开整个文本文件的情况下将星号添加到行中.鉴于我们必须重复该过程未知次数,这可能非常有用.

我想测试一堆形式的基因组位置:

chr4:154723876-154724615
chr6:139580853-139581090
chr18:30440532-30441569

我想看看它们是位于UTR还是内含子或外显子或基因间序列.我不关心这些坐标是哪些基因的内含子(等)的信息.

我假设每个已知的遗传元件(如外显子)都定义了基因组位置(每条染色体上基因组的起始位置).我知道外显子和内含子也是如此,例如Ensembl在基因组中有每个外显子的ID:参见Mus musclulus中Amy1基因的外显子和内含子的例子.我想用上面的位置列表查询这些位置的数据库,如果两者之间有重叠(理想情况下我应该能够指定重叠,比如说,至少10bp,但如果不是,我可以) ,我应该受欢迎(是的,这个区域在外显子/内含子/)

差点在于我有几千个这样的位置,并且理想情况下想要一次性查询它们并且作为输出有一个表格,其中每个位置将被分配为"内含子/外显子/ utr/intergenic".有机体是Mus musculus,位置来自整个基因组.

我现在不能提供我正在尝试做的代码示例,因为我不知道从哪里开始 - 如果我有一个包或任何内容可以帮助我找到解决方案.

如果我可以在R中完成,那将是完美的,但AFAIK我不能在biomaRt中做到这一点,我找不到一个包来做它.我想到了Galaxy,但是考虑到他们这样做的非平凡方式和他们产生的奇怪输出,我宁愿坚持R.你知道的魔鬼等.

非常感谢帮助.

还是一个Python新手所以请放轻松我...

我有一个字典设置:

new_dict

我想过滤以返回键,其中每个键附加的任何值与我设置的现有列表中的值匹配:

list(data.Mapped_gene)

有任何想法吗？

编辑:我仍然无法完成这项工作.

如果有帮助,csv表和键都是字符串.

以下是扩大理解的完整代码:

import csv    
new_dict = {}
with open(raw_input("Enter csv file (including path)"), 'rb') as f:
  reader = csv.reader(f)
  for row in reader:
    if row[0] in new_dict:
      new_dict[row[0]].append(row[1:])
    else:
      new_dict[row[0]] = row[1:]
print new_dict

#modified from: http://bit.ly/1iOS7Gu
import pandas
colnames = ['Date Added to Catalog',    'PUBMEDID', 'First Author', 'Date',     'Journal',  'Link', 'Study',    'DT',   'Initial Sample Size',  'Replication Sample Size',  'Region',   'Chr_id',   'Chr_pos',  'Reported Gene(s)', 'Mapped_gene',  'p-Value',  'Pvalue_mlog',  'p-Value (text)',   'OR or beta',   '95% …

我有这个VCF格式的文件，我想在R中读取此文件。但是，此文件包含一些我想跳过的多余行。我想在行以匹配行开始的结果中得到类似的结果#CHROM。

这是我尝试过的：

chromo1<-try(scan(myfile.vcf,what=character(),n=5000,sep="\n",skip=0,fill=TRUE,na.strings="",quote="\"")) ## find the start of the vcf file
skip.lines<-grep("^#CHROM",chromo1)


column.labels<-read.delim(myfile.vcf,header=F,nrows=1,skip=(skip.lines-1),sep="\t",fill=TRUE,stringsAsFactors=FALSE,na.strings="",quote="\"")
num.vars<-dim(column.labels)[2]

myfile.vcf

    #not wanted line
    #unnecessary line
    #junk line
    #CHROM  POS     ID      REF     ALT
    11      33443   3        A       T
    12      33445   5        A       G

结果

    #CHROM  POS     ID      REF     ALT
    11      33443   3        A       T
    12      33445   5        A       G

这可能是一个非常奇怪的问题,它肯定是.我不太熟悉编程语言是如何用传统方法制作的,所以我想知道,是否可以设计无语法编程语言？这意味着任何输入都是有效的并执行某个计算,并且相同的输入将始终执行相同的操作.将不会出现语法错误(允许逻辑和运行时错误,程序可能崩溃,进行随机计算等).

我想到了这一点,因为遗传学基本上是我的理解,就像那样.

编辑:我认为存在一些误解.无语法只是意味着所有输入都将计算,解释器/编译程序将遵循该特定指令集,但可能是随机的.

此外,它必须匹配每个输入都有1个且只有1个输出的事实.诸如语法错误之类的内容违反了该规则.

编辑2很多人都在使用语法部分.忘记语法,关注任何输入将产生UNIQUE输出的事实.

snp1 <- c("AA", "AT", "AA", "TT", "AA", "AT", "AA", "AA", "AA", "AT")
snp2 <- c("GG", "GC", "GG", "CC", "CC", "GC", "GG", "GG", "GG", "GC")
df1 <- data.frame(snp1, snp2)

num1 <- c(1, 2, 1, 3, 1, 2, 1, 1, 1, 2)
num2 <- c(1, 2, 1, 3, 3, 2, 1, 1, 1, 2)
df2 <- data.frame(num1, num2)

这是在R中完成的.我有一个对象df1,我想将其转换为df2.对于df1中的每一列,最常见的值转换为1,第二个最常见的值转换为2,等等.我该如何有效地做到这一点？

我正在尝试在 R 中转置数据帧，但运气不佳。

数据框包含一个表观遗传数据集，第一列中有 300,000 多个 CpG 位点。74 个额外的列被分成对照组和实验组（癌症 = 69，正常 = 5）。

我已将数据框转换为矩阵，以便可以转置数据并将其转换回数值。但是，每次我尝试将数据转换回数据框时，它最终都会变成一个列表。

我正在尝试执行以下操作：

1. 将第一列设为标题。

2. 将标题的名称设为第一列中的值（我想删除的 X 除外）。

   数据[1:10, 1:10]

       X Tumor.33 Normal.01 Tumor.01 Tumor.34 Tumor.35 Tumor.02 Tumor.03
   
   
   cg00000029  0.29224   0.32605  0.58762  0.32397  0.23482  0.24012  0.22941
   
   cg00000108  0.91243   0.89785  0.92337  0.90080  0.91220  0.92256  0.92709
   
   cg00000109  0.77676   0.73910  0.81545  0.73603  0.76276  0.85808  0.85142
   
   cg00000165  0.34261   0.30960  0.56392  0.32363  0.33980  0.61755  0.70855
   
   cg00000236  0.84688   0.80654  0.84423  0.80935  0.85600  0.87766  0.83509
   
   cg00000289  0.61535   0.60874  0.62496  0.66421  0.60556  0.66824  0.65243
   
   cg00000292  0.72491   0.63333 …

   Run Code Online (Sandbox Code Playgroud)

我正在尝试使用{pegas}的haploNet功能来绘制单倍型网络,但是我在同一个饼图中将来自不同群体的相同单倍型变得困难.我可以用以下脚本构建单倍型网:

x <- read.dna(file="x.fas",format="fasta")
h <- haplotype(x)
net <- haploNet(h)
plot(net)

我想在dnabin数据中设置每个分类群的原始种群的标签,因此我可以在得到的网络中得到不同颜色的饼图(来自不同种群的单倍型).我还想删除生成的单倍型网络中的重叠圆圈.

谢谢你的帮助!

一个例子:

> data(woodmouse)
> x <- woodmouse[sample(15, size = 110, replace = TRUE), ]
> h <- haplotype(x)
> net <- haploNet(h)
> plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.8)

该脚本用于使用{pegas}构建单倍型网络.较大的圆圈代表某种类型的更多单倍型.我想知道如何在dnabin矩阵中设置单倍型的起源,因此它们会在网络中以不同的颜色出现.

当使用 haploNet 包在单倍型网络上绘制一些图时，我使用了 Internet 上可用的脚本来执行此操作。不过我觉得有什么地方不对。该脚本以 woodmouse 示例的形式提供。我使用的代码是：

x <- read.dna(file="Masto.fasta",format="fasta")
h <- haplotype(x)
net <- haploNet(h)
plot(net)

plot(net, size = attr(net, "freq"), fast = TRUE)
plot(net, size = attr(net, "freq"))
plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.8

table(rownames(x))

ind.hap<-with(
    stack(setNames(attr(h, "index"), rownames(h))), 
    table(hap=ind, pop=rownames(x)[values])
)
ind.hap 

plot(net, size=attr(net, "freq"), scale.ratio = 2, cex = 0.8, pie=ind.hap)
legend(50,50, colnames(ind.hap), col=rainbow(ncol(ind.hap)), pch=20)

legend(x=7,y=10,c("Baeti ero","Felege weyni","Golgole naele","Hagare selam","Ruba feleg","Ziway"),c("red","yellow","green","turquoise","blue","magenta"))

但是，在绘制 ind.hap 时，您会注意到某些行不在正确的位置。你可以在这里看到这个：

      pop
hap    Baetiero ETH022 ETH742 Felegeweyni Golgolenaele Rubafeleg
  I …

在遗传学中,非常小的p值是常见的(例如10 ^ -400),当我在R中得分很大时,我正在寻找一种获得非常小的p值(双尾)的方法,例如:

z=40
pvalue = 2*pnorm(abs(z), lower.tail = F)

这给了我零而不是非常小的值,这是非常重要的.

我想使用R和sub提取符号周围的字符.我尝试了许多正则表达式,但我没有得到我想要的东西.

我的矢量:

c("G>GA", "T>A", "G>A", "G>A", "A>T", "CT>C", "T>C", "T>C", "A>T", "T>C", "T>A", "A>G", "CCGCCGCGGCCGCCGTCTTCCACCAACAACATGGCGGA>C", "C>T", "T>A", "T>C", "T>G", "G>C", "T>G", "T>A", "G>A")

我之前和之后只需要一个角色>.

我最好的尝试是:

sub("(.*?)>", ">", aa, perl = TRUE)

我有一个 .fam、.bed 和 .bim 文件，其中包含少数个人的标记。我需要将其转换为 VCF 文件。

有人可以帮忙创建一个 VCF 文件吗？有没有开源工具可以做到这一点？

我正在寻找一种方法来查找 python 中两个字符串之间不匹配的总数。我的输入是一个看起来像这样的列表

['sequence=AGATGG', 'sequence=AGCTAG', 'sequence=TGCTAG',
 'sequence=AGGTAG', 'sequence=AGCTAG', 'sequence=AGAGAG']

对于每个字符串，我想看看它与序列有多少差异"sequence=AGATAA"。因此，如果输入[0]来自上面的列表，则输出将如下所示：

sequence=AGATGG, 2

我不知道是否将每个字母拆分为单独的列表，或者是否应该尝试以某种方式比较整个字符串。任何帮助都是有用的，谢谢