标签: bioinformatics

我想做一个相对简单的情节(让人想起如下的时间线:http://www.ats.ucla.edu/stat/sas/code/timeline.gif),但不是时间在x轴上,它将成为基因组中的基础位置."时间跨度"将是DNA序列支架的覆盖距离,显示它们落入基因组的范围,它们重叠的位置和没有覆盖的位置.这是我正在寻找的粗略模型,显示rRNA的重叠群覆盖,(我遗漏了,但需要,x轴显示开始和停止的位置,以及重叠群的标记(彩色线)): http://i.imgur.com/MDABx.png,具有以下坐标:

Contig# Start1  Stop1   Start2  Stop2   Start3  Stop3   Start4  Stop4
1   1   90  90  100 120 150 200 400
2   1   100 120 150 200 400 NA  NA
3   1   30  90  100 120 135 200 400
4   1   100 120 140 200 400 NA  NA
5   -35 80  90  100 130 150 200 400
6   1   100 200 300 360 400 NA  NA

我很确定这可以在R中完成,可能使用ggplot2,但由于某种原因我无法弄明白.

我有一张桌子,开头就在这里:

>df
TargetID    SM_H1455_121005_4   SM_H1456_121005_1   SM_H1457_121005_7   SM_H1461_121005_2   SM_H1462_121005_8   SM_H1463_121005_1   SM_K1566_121005_6   

ENSG00000002549.7   2286    2468    2498    1696    2044    11536   4100    5460        
ENSG00000002587.5   10  0   6   2   0   2   34  
ENSG00000002726.15  8   14  0   2   16  2   4   
ENSG00000002745.8   6   2   2   0   0   4   6   

我想浏览每一行并删除包含0的每一行.这是我的代码:

for(i in 2:length(df)) {
df2 <- df[df[,i]!=0,]
}

但是这段代码不会删除其中包含0的所有行.它删除了一些行,但我不确定为什么,我不认为这是正确的.有什么建议？

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

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

我想将多个指定位置的值插入到 python 字符串/数组中。

例如我的输入字符串：SARLSAMLVPVTPEVKPK

在指定位置：1,5,12

所需的输出：S*ARLS*AMLVPVT*PEVKPK

我试过：

seq="SARLSAMLVPVTPEVKPK" #string
pos=[1,5,12] #positions
arr=list(seq) #convert string to array
arr.insert(pos,"*") # NOT WORK!
arr.insert(pos[0],"*")
print(''.join(arr))

看来我一次只能插入一个位置，因此下一次插入的指定位置的索引必须更改。有没有一种优雅的方法来做到这一点，或者我是否必须循环遍历插入位置，为每个附加插入位置添加+1？我希望这是有道理的！

非常感谢，卷毛。

我有 DNA 序列数据。例如，

X="ACGGGT"
Y="ACGGT"

我想知道对齐分数，因此我使用了biopythonpairwise2函数。例如，

from Bio import pairwise2
from Bio.pairwise2 import format_alignment

alignments = pairwise2.align.globalxx(X, Y)
for a in alignments:
    print(format_alignment(*a))

这成功地显示了 DNA 比对，但我只需要如下的分数。有没有办法只显示分数？

我使用了biopython，但如果有更好的方法，我们将不胜感激。

我希望提高我的生物信息学算法的速度，这需要比较“A”、“C”、“G”、“T”之一的字符（例如计算“A”==“C”）

由于字符的大小为 8 位，因此在最坏的情况下需要对二进制数进行 8 次比较。我的猜测是，通过将 'A'、'C'、'G'、'T' 表示为一对二进制数（例如，将 'A' 表示为 make_pair(false,false)），我想我可以提高速度3~4 次，因为它现在最多只需要 2 次二进制比较。

我尝试使用一对布尔值，但速度实际上下降了大约 30%。

表示四个字符和计算相等性的最快方法是什么？内存使用对我来说并不是什么大问题。

供您参考，我使用的是 C++11 编译器。先感谢您。

我有一个 ID 数据集：

VARIANT_ID
01_1254436_A_G_1
02_2254436_A_G_1 
03_3255436_A_G_1 
10_10344745_A_G_1 
11_11256437_A_G_1 
11_11343426_A_G_1 
12_12222431_A_G_1
14_14200436_A_G_1 
15_15256789_A_G_1 

我希望仅从 ID 以 01-09 开头的行的开头字符中删除 0，但在不进一步删除列中的其他 0 的情况下，我很难做到这一点，并且只能看到其他语言的类似问题。我想要的输出是：

VARIANT_ID
1_1254436_A_G_1
2_2254436_A_G_1 
3_3255436_A_G_1 
10_10344745_A_G_1 
11_11256437_A_G_1 
11_11343426_A_G_1 
12_12222431_A_G_1
14_14200436_A_G_1 
15_15256789_A_G_1 

仅删除了每行开头的零，如何指定这一点？我来自生物学背景，因此任何帮助将不胜感激。

输入数据：

structure(list(VARIANT_ID = c("01_1254436_A_G_1", "02_2254436_A_G_1", 
"03_3255436_A_G_1", "10_10344745_A_G_1", "11_11256437_A_G_1", 
"11_11343426_A_G_1", "12_12222431_A_G_1", "14_14200436_A_G_1", 
"15_15256789_A_G_1")), row.names = c(NA, -9L), class = c("data.table", 
"data.frame"))

我正在尝试从 GEO/SRA 下载 BAM 格式的数据集，我可以使用它在 RStudio 中进行分析。

我尝试使用这种方法：我下载了 .sra 并将其转换为 .bam

prefetch GSM269238
sam-dump C:\Users\Desktop\sratoolkit.2.10.8-win64\bin\ncbi\SRA\sra\GSM2692389.sra --output-file GSM2692389.bam

然而，在 RStudio 中这不起作用，并返回一个错误，说它无法读取 bam 文件这是我的 R 代码；我正在使用 RSamTools

> bamfiles <- list.files("directory redacted due to privacy", ".bam")
> file.exists(bamfiles)
[1] TRUE
> 
> 
> #---> Define bam files for count step on Rsamtools
> 
> library("Rsamtools")
> bamfiles <- BamFileList(bamfiles, yieldSize=2000000)
> seqinfo(bamfiles)
Error in value[[3L]](cond) : 
  failed to open BamFile: SAM/BAM header missing or empty
  file: 'GSM2692389.bam'

有谁知道如何帮助我将 SRA 数据下载到可读的 .bam 文件中？任何帮助或指导将不胜感激，因为我真的很努力在截止日期前完成这件事。

我正在使用cellranger mkref并面临 GTF（自定义 gtf 文件）的奇怪 python 问题：

Traceback (most recent call last):
  File "/home/user/cellranger-6.0.1/lib/python/cellranger/reference.py", line 750, in validate_gtf
    subprocess.check_output(cmd, stderr=subprocess.STDOUT)
  File "/home/user/cellranger-6.0.1/external/anaconda/lib/python3.7/subprocess.py", line 411, in check_output
    **kwargs).stdout
  File "/home/user/cellranger-6.0.1/external/anaconda/lib/python3.7/subprocess.py", line 512, in run
    output=stdout, stderr=stderr)
subprocess.CalledProcessError: Command '['gtf_to_gene_index', '/home/user/cellranger-6.0.1/indexes', '/home/user/cellranger-6.0.1/indexes/tmp74f_vsxg.json']' returned non-zero exit status 1.

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "/home/user/cellranger-6.0.1/bin/rna/mkref", line 139, in <module>
    main()
  File "/home/user/cellranger-6.0.1/bin/rna/mkref", line 130, in main
    reference_builder.build_gex_reference()
  File "/home/user/cellranger-6.0.1/lib/python/cellranger/reference.py", line 613, in …

我对我的问题有一个可行的解决方案，但速度很慢。我很好奇推荐的加速方法，并想看看它能达到多快。这是一个示例输入文件

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   0|1:7,28:35:99:0|1:616504_T_C:787,0,177:616504   
Chr16 616504 T     C     BESC.1    0/0:48,0:48:99:.:.:0,114,1710:.                  
Chr16 616504 T     C     BESC.10   1|1:0,23:23:72:1|1:616504_T_C:1059,72,0:616504   
Chr16 616504 T     C     BESC.100  0/0:34,0:34:96:.:.:0,96,1440:.                   
Chr16 616504 T     C     BESC.1001 0/0:47,0:47:99:.:.:0,120,1800:.                  
Chr16 616504 T     C     BESC.1002 0/0:39,0:39:99:.:.:0,108,948:.    

           

目标是从value列中取出第一个和第三个字符并对它们求和，然后输出一个类似的文件，其中值列替换为该总和。前两行的示例输出：

CHROM    POS REF   ALT   Geno      value                                                                                       
Chr16 616504 T     C     X93.968   1   
Chr16 616504 T     C     BESC.1    0   

这是我当前的解决方案，其中 STDIN 1 是输入文件名，STDIN 2 是输出文件名：

#!/bin/bash
i=0
len=$(cat $1 | wc -l)

touch $2
while …