标签: bioinformatics

我正在研究一个python项目,在那里我研究RNA结构演化(例如表示为字符串:"(((...)))"其中括号代表碱基对.关键是我有一个理想的结构和一个向理想结构发展的人口.我实现了一切,但是我想添加一个功能,我可以获得"桶数",即每代人口中最具代表性的k结构.

我在考虑使用k-means算法,但我不确定如何将它与字符串一起使用.我找到了scipy.cluster.vq,但我不知道如何在我的情况下使用它.

谢谢!

哪些函数式编程语言可以轻松获得生物信息库？

(不要包含Ruby等多范式语言)

更新:列出哪些主要的函数式编程语言目前不能轻松访问生物信息库也是受欢迎的.

我有一个小的DNA序列快速文件,看起来像这样:

>NM_000016 700 200 234
ACATATTGGAGGCCGAAACAATGAGGCGTGATCAACTCAGTATATCAC

>NM_000775 700 124 236
CTAACCTCTCCCAGTGTGGAACCTCTATCTCATGAGAAAGCTGGGATGAG

>NM_003820 700 111 222
ATTTCCTCCTGCTGCCCGGGAGGTAACACCCTGGACCCCTGGAGTCTGCA

问题:

1)如何将此fasta文件读入R作为数据帧,其中每一行是序列记录,第一列是refseqID,第二列是序列.

2)如何在(开始,结束)位置提取子序列？

NM_000016 1  3 #"ACA"
NM_000775 2  6 #"TAACC"
NM_003820 3  5 #"TTC"

我正在使用python/pysam来分析测序数据.在其教程(pysam - 读取和编写SAM文件的界面)中,命令配合说:

'这种方法对于高吞吐量处理来说太慢了.如果需要使用其配合处理读取,则从读取名称的排序文件开始工作,或者更好的是缓存读取.

你会如何"缓存读取"？

我正在尝试解决Rosalind在给定序列中计算核苷酸的基本问题,并将结果返回到列表中.对于那些不熟悉生物信息学的人来说,它只计算一个字符串中4个不同字符('A','C','G','T')的出现次数.

我希望collections.Counter这是最快的方法(首先是因为他们声称是高性能,第二是因为我看到很多人使用它来解决这个特定的问题).

但令我惊讶的是这种方法最慢!

我比较了三种不同的方法,使用timeit和运行两种类型的实验:

• 长时间运行几次
• 多次运行短序列.

这是我的代码:

import timeit
from collections import Counter

# Method1: using count
def method1(seq):
    return [seq.count('A'), seq.count('C'), seq.count('G'), seq.count('T')]

# method 2: using a loop
def method2(seq):
    r = [0, 0, 0, 0]
    for i in seq:
        if i == 'A':
            r[0] += 1
        elif i == 'C':
            r[1] += 1
        elif i == 'G':
            r[2] += 1
        else:
            r[3] += 1
    return r

# method 3: using Collections.counter …

我正在寻找一种更快的方法来计算从FASTA文件读入的DNA字符串的GC内容.这归结为取一个字符串并计算字母'G'或'C'出现的次数.我还想指定要考虑的字符范围.

我有一个相当慢的工作函数,它导致我的代码瓶颈.它看起来像这样:

##
## count the number of GCs in the characters between start and stop
##
gcCount <-  function(line, st, sp){
  chars = strsplit(as.character(line),"")[[1]]
  numGC = 0
  for(j in st:sp){
    ##nested ifs faster than an OR (|) construction
    if(chars[[j]] == "g"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "G"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "c"){
      numGC <- numGC + 1
    }else if(chars[[j]] == "C"){
      numGC <- numGC + 1
    }
  }
  return(numGC)
}

运行Rprof给我以下输出:

> …

我有两个data.frames,每个都有三列:chrom,start和stop,让我们称它们为rangesA和rangesB.对于每行rangeA,我想找到rangeB中哪一行(如果有的话)完全包含rangesA行 - 我的意思是rangesAChrom == rangesBChrom, rangesAStart >= rangesBStart and rangesAStop <= rangesBStop.

现在我正在做以下,我不太喜欢.请注意,由于其他原因,我正在循环遍历rangeA的行,但这些原因都不是什么大问题,它最终会让这些特定解决方案更具可读性.

rangesA:

chrom   start   stop
 5       100     105
 1       200     250
 9       275     300

rangesB:

chrom    start    stop
  1       200      265
  5       99       106
  9       275      290

对于范围A中的每一行:

matches <- which((rangesB[,'chrom']  == rangesA[row,'chrom']) &&
                 (rangesB[,'start'] <= rangesA[row, 'start']) &&
                 (rangesB[,'stop'] >= rangesA[row, 'stop']))

我认为必须有一个更好的(并且更好,我的意思是比大范围A和rangeB的大型实例更快)这样做的方式比循环这个结构.有任何想法吗？

我几乎让我的needleman-wunsch实现工作,但我对如何处理特定情况下的追溯感到困惑.

我们的想法是,为了重新构建序列(最长路径),我们重新计算以确定得分来自的矩阵.我遇到问题的边缘情况是右下方得分不在匹配矩阵中,但是在插入列矩阵中(意味着得到的跟踪序列应该有插入).

这些序列以a2m格式记录,其中序列中的插入被记录为小写字符.所以在最终输出中ZZ,AAAC应该是对齐AAac.当我用手追溯时,我最终会AAAc因为我只访问Ic矩阵一次.这是我的白板图片.如你所见,我有三个黑色箭头和一个绿色箭头,这就是我追溯给我的原因AAAc.我应该计算第一个单元格,然后停在1,1位置？我不确定如何改变我实现这一点的方式.

注意,这里使用的替换矩阵是BLOSUM62.复发关系是

M(i,j) = max(Ic(i-1,j-1)+subst, M(i-1,j-1)+subst, Ir(i-1,j-1)+subst)
Ic(i,j) = max(Ic(i,j-1)-extend, M(i,j-1)-open, Ir(i,j-1)-double)
Ir(i,j) = max(Ic(i-1,j)-double, M(i-1,j)-open, Ir(i-1,j)-extend)

编辑:这是traceback_col_seq函数重写为更清洁.请注意,score_cell现在返回thisM,thisC,thisR而不是其中的最大值.这个版本将对齐评为AaAc,仍然存在同样的问题,现在又出现了另一个问题,即为什么它会在1,2处再次进入Ic.但是这个版本更清晰.

def traceback_col_seq(self):
    i, j = self.maxI-1, self.maxJ-1
    self.traceback = list()
    matrixDict = {0:'M',1:'Ic',2:'Ir',3:'M',4:'Ic',5:'Ir',6:'M',7:'Ic',8:'Ir'}
    while i > 0 or j > 0:
        chars = self.col_seq[j-1] + self.row_seq[i-1]
        thisM, thisC, thisR = self.score_cell(i, j, chars)
        cell = thisM + thisC + thisR
        prevMatrix = matrixDict[cell.index(max(cell))]
        print(cell, prevMatrix,i,j)
        if prevMatrix == …

我有一个Newick树是通过比较4-9bp长DNA序列的推定DNA调节基序的位置权重矩阵(PWM或PSSM)的相似性(欧几里德距离)而构建的.

这个树的交互式版本在iTol(这里)上,您可以自由使用 - 只需在设置参数后按"更新树":

我的具体目标是:如果它们与最近的父进化枝的平均距离是<X(ETE2 Python包),则将这些图案(提示/终端节点/叶子)折叠在一起.这在生物学上是有意义的,因为一些基因调控DNA基序可以彼此同源(旁系同源物或直向同源物).这种折叠可以通过上面链接的iTol GUI完成,例如,如果你选择X = 0.001,那么一些图案会折叠成三角形(图案族).

我的问题:任何人都可以提出一种算法,可以输出或帮助可视化X的哪个值适合"最大化折叠图案的生物或统计相关性"？理想情况下,当对X绘制时,树的某些属性会有一些明显的阶跃变化,这表明算法是一个合理的X.是否有任何已知的算法/脚本/包？也许代码会针对X的值绘制一些统计数据？我已经尝试绘制X与平均簇大小(matplotlib),但我没有看到明显的"步骤增加"来通知我使用X的值:

我的代码和数据:我的Python脚本的链接是[这里] [8],我已经对它进行了大量评论,它将为您生成树数据和上图(使用参数d_from,d_to和d_step来探索距离切割-offs,X).如果你有easy-install和Python,你只需执行这两个bash命令就需要安装ete2:

apt-get install python-setuptools python-numpy python-qt4 python-scipy python-mysqldb python-lxml

easy_install -U ete2

我目前正在使用字符串形式的DNA序列,其中内含子是小写字母,外显子是大写字母.该方法的目的是尽可能快地以String的形式检索外显子.

序列示例:

ATGGATGACAGgtgagaggacactcgggtcccagccccaggctctgccctcaggaagggggtcagctctcaggggcatctccctctcacagcccagccctggggatgatgtgggagccccatttatacacggtgcctccttctctcctagAGCCTACATAG

我的第一个版本是使用String replaceAll()方法,但它特别慢:

public String getExons(String sequence) {
    return sequence.replaceAll("[atgc]", "");
}

所以我尝试了一个新版本,它提高了性能,但仍然很慢:

public String getExons(String sequence) {
    StringBuilder exonBuilder = new StringBuilder();

    for (int i = 0; i < sequence.length(); i++) {
        char c = sequence.charAt(i);
        if (c == 'A' || c == 'T' || c == 'G' || c == 'C') exonBuilder.append(c);
    }
    return exonBuilder.toString();

还有另一种方法可以提高性能吗？