标签: bioinformatics

我有一个fasta文件,如下所示.我想将三个字母的代码转换为一个字母代码.我怎么能用python或R做到这一点？

>2ppo
ARGHISLEULEULYS
>3oot
METHISARGARGMET

期望的输出

>2ppo
RHLLK
>3oot
MHRRM

你的建议将不胜感激!

我的数据按V6中的ID分组,按位置排序(V1:V3):

dt
      V1      V2      V3 V4 V5                 V6
 1: chr1 3054233 3054733  .  + ENSMUSG00000090025
 2: chr1 3102016 3102125  .  + ENSMUSG00000064842
 3: chr1 3205901 3207317  .  - ENSMUSG00000051951
 4: chr1 3206523 3207317  .  - ENSMUSG00000051951
 5: chr1 3213439 3215632  .  - ENSMUSG00000051951
 6: chr1 3213609 3216344  .  - ENSMUSG00000051951
 7: chr1 3214482 3216968  .  - ENSMUSG00000051951
 8: chr1 3421702 3421901  .  - ENSMUSG00000051951
 9: chr1 3466587 3466687  .  + ENSMUSG00000089699
10: chr1 3513405 3513553  .  + ENSMUSG00000089699

我想要做的是添加一个带位置索引的额外列,也就是说,每个组在V6中第一个元素是"1",第二个元素是"2",依此类推.我可以使用ddply和自定义函数实现:

rankExons …

我有一个包含Ensembl ID的data.frame在一列中; 我想为该列的值找到相应的基因符号,并将它们添加到我的数据框中的新列.我使用bioMaRt但它找不到任何Ensembl ID!

这是我的示例数据(df[1:2,]):

row.names organism    gene
41  Homo-Sapiens ENSP00000335357
115 Homo-Sapiens ENSP00000227378

我希望得到这样的东西

row.names organism    gene         id
41  Homo-Sapiens ENSP00000335357   CDKN3
115 Homo-Sapiens ENSP00000227378   HSPA8

这是我的代码:

library('biomaRt')
mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))
genes <- df$genes
df$id <- NA
G_list <- getBM(filters= "ensembl_gene_id", attributes= c("ensembl_gene_id",
"entrezgene", "description"),values=genes,mart= mart)

然后,当我检查G_list时,我得到了这个

[1] ensembl_gene_id entrezgene      description  <0 rows> (or 0-length row.names)

所以我无法将G_list添加到我的df中!因为没有什么可补充的!

提前致谢,

参与生物信息学的人应该知道哪些数据结构？我想,任何人都应该知道的名单,哈希,平衡树等,但我希望有特定领域的数据结构.有没有专门讨论这个主题的书？

我试图找到一种方法来折叠具有相交范围的行,用"开始"和"停止"列表示,并将折叠值记录到新列中.例如,我有这个数据框:

my.df<- data.frame(chrom=c(1,1,1,1,14,16,16), name=c("a","b","c","d","e","f","g"), start=as.numeric(c(0,70001,70203,70060, 40004, 50000872, 50000872)), stop=as.numeric(c(71200,71200,80001,71051, 42004, 50000890, 51000952)))


chrom name  start   stop
 1    a        0    71200
 1    b    70001    71200
 1    c    70203    80001
 1    d    70060    71051
14    e    40004    42004
16    f 50000872 50000890
16    g 50000872 51000952

我试图找到重叠的范围并记录"开始"和"停止"中折叠的重叠行所涵盖的最大范围以及折叠行的名称,所以我会得到:

chrom start   stop      name
 1    70001    80001    a,b,c,d
14    40004    42004    e
16    50000872 51000952 f,g

我想我可以像这样使用IRanges包:

library(IRanges)
ranges <- split(IRanges(my.df$start, my.df$stop), my.df$chrom)

但是我在收到折叠列时遇到了麻烦:我尝试过findOvarlaps但是这个

ov <- findOverlaps(ranges, ranges, type="any")

但我不认为这是对的.

任何帮助将非常感激.

编码新手.Pytho/biopython新手; 这是我在网上的第一个问题.如何打开压缩的fasta.gz文件以提取信息并在我的函数中执行计算.这是我正在尝试做的一个简化示例(我尝试过不同的方法),以及错误是什么.我正在使用的gzip命令似乎不起作用.

with gzip.open("practicezip.fasta.gz", "r") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

Traceback (most recent call last):

  File "<ipython-input-192-a94ad3309a16>", line 2, in <module>
    for record in SeqIO.parse(handle, "fasta"):

  File "C:\Users\Anaconda3\lib\site-packages\Bio\SeqIO\__init__.py", line 600, in parse
    for r in i:

  File "C:\Users\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py", line 122, in FastaIterator
    for title, sequence in SimpleFastaParser(handle):

  File "C:\Users\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py", line 46, in SimpleFastaParser
    if line[0] == ">":

IndexError: index out of range

您正在使用基因表达矩阵,片段计数来计算差异表达的基因.我想知道如何删除值为0的行.然后我的数据集将是紧凑的,并且我将使用此矩阵为下游分析提供更少的虚假结果.

输入

gene    ZPT.1   ZPT.0   ZPT.2   ZPT.3   PDGT.1  PDGT.0
XLOC_000001 3516    626 1277    770 4309    9030
XLOC_000002 342 82  185 72  835 1095
XLOC_000003 2000    361 867 438 454 687
XLOC_000004 143 30  67  37  90  236
XLOC_000005 0   0   0   0   0   0
XLOC_000006 0   0   0   0   0   0
XLOC_000007 0   0   0   0   1   3
XLOC_000008 0   0   0   0   0   0
XLOC_000009 0   0   0   0   0   0
XLOC_000010 7   1   5   3   0   1
XLOC_000011 …

我正在尝试使用仿射空位罚函数实现用于局部序列比对的Smith-Waterman算法.我想我理解如何启动和计算计算校准分数所需的矩阵,但我不知道如何追溯以找到对齐.要生成所需的3个矩阵,我有以下代码

for j in range(1, len2):
    for i in range(1, len1):
        fxOpen = F[i][j-1] + gap
        xExtend = Ix[i][j-1] + extend
        Ix[i][j] = max(fxOpen, xExtend)

        fyOpen = F[i-1][j] + gap
        yExtend = Iy[i-1][j] + extend
        Iy[i][j] = max(fyOpen, yExtend)

        matchScore = (F[i-1][j-1]  + simMatrixDict[seq1[i-1]+seq2[j-1]])
        xScore = Ix[i-1][j-1] + simMatrixDict[seq1[i-1]+seq2[j-1]]
        yScore = Iy[i-1][j-1] + simMatrixDict[seq1[i-1]+seq2[j-1]]
        F[i][j] = max(0, matchScore, xScore, yScore)

我不确定我是否需要一个矩阵用于追溯,或者只需要1？任何关于如何从F中的最大分数追溯的澄清将非常感激.

我有几个SNP ID(即rs16828074,rs17232800等),我想在UCSC基因组网站的Hg19基因组中找到它们的坐标.

我更愿意用它R来实现这个目标.怎么做？

我有一个Fortran程序,并希望在python中为多个文件执行它.我有2000个输入文件,但在我的Fortran代码中,我一次只能运行一个文件.我应该如何在python中调用Fortran程序？

我的剧本:

import subprocess
import glob
input = glob.glob('C:/Users/Vishnu/Desktop/Fortran_Program_Rum/*.txt')
output = glob.glob('C:/Users/Vishnu/Desktop/Fortran_Program_Rum/Output/')
f = open("output", "w")
for i in input:
    subprocess.Popen(["FORTRAN ~/C:/Users/Vishnu/Desktop/Fortran_Program_Rum/phase1.f", "--domain "+i])
f.write(i)

错误:

runfile('C:/Users/Vishnu/Desktop/test_fn/test.py', wdir='C:/Users/Vishnu/Desktop/test_fn')
Traceback (most recent call last):

   File "<ipython-input-3-f8f378816004>", line 1, in <module>
runfile('C:/Users/Vishnu/Desktop/test_fn/test.py', wdir='C:/Users/Vishnu/Desktop/test_fn')

  File "C:\Users\Vishnu\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 866, in runfile
execfile(filename, namespace)

  File "C:\Users\Vishnu\Anaconda3\lib\site-packages\spyder\utils\site\sitecustomize.py", line 102, in execfile
exec(compile(f.read(), filename, 'exec'), namespace)

  File "C:/Users/Vishnu/Desktop/test_fn/test.py", line 30, in <module>
subprocess.Popen(["FORTRAN ~/C:/Users/Vishnu/Desktop/Fortran_Program_Rum/phase1.f", "--domain "+i])

  File "C:\Users\Vishnu\Anaconda3\lib\subprocess.py", line 707, in __init__
restore_signals, start_new_session)

  File "C:\Users\Vishnu\Anaconda3\lib\subprocess.py", line …