我有大约 500 个 fasta 格式的蛋白质序列，我是从 Blastp 搜索中得到的。从这些序列中，我需要有蛋白质名称、生物体、Uniprot ID 以及蛋白质家族（如果可能），以便我可以用这些信息构建一个表格。
有什么办法可以用python来做吗？一些与 Uniprot 通信的功能？如何解析 fasta 标头中的信息？
我是 R 初学者。我正在使用 R 来分析我的大型下一代测序 vcf 文件，但遇到了一些困难。我已将非常大的 vcf 文件作为数据框（2446824 个观测值，共 177 个变量）导入，并仅使用我感兴趣的 3 个样本（2446824 个观测值，共 29 个变量）制作了一个子集。
我现在希望进一步减少尺寸（将行数减少到 200000 左右）。我一直在尝试使用 grep，但无法让它工作。我得到的错误是
Error in "0/1" | "1/0" :
operations are possible only for numeric, logical or complex types
这是我正在使用的文件的一个小示例部分。
Chr Start End Ref Alt Func.refGene INFO FORMAT Run.Sample1 Run.Sample2 Run.Sample3
489 1 909221 909221 T C PASS GT:AD:DP:GQ:PL 0/1:11,0:11:33:0,33,381 ./. ./.
490 1 909238 909238 G C PASS GT:AD:DP:GQ:PL 0/1:11,6:17:99:171,0,274 0/1:6,5:11:99:159,0,116 1/1:0,15:15:36:441,36,0
491 1 909242 909242 A G PASS GT:AD:DP:GQ:PL 0/1:16,4:13:45:0,45,532 …有人知道如何仅使用GenBank代码登录名和biopython从GenBank中的数据中获取科学名称（或所有功能）。例如：
>>> From Bio import Entrez
>>> Entrez.email = someuser@mail.com
>>> Input = Entrez.someFunction(db="nucleotide", term="AY851612")
>>> output = Entrez.read(Input)
>>> print output
"Austrocylindropuntia subulata"
还是：
>>> print output
"LOCUS AY851612 892 bp DNA linear PLN 10-APR-2007
DEFINITION Opuntia subulata rpl16 gene, intron; chloroplast.
ACCESSION AY851612
VERSION AY851612.1 GI:57240072
KEYWORDS .
SOURCE chloroplast Austrocylindropuntia subulata
ORGANISM Austrocylindropuntia subulata
Eukaryota; Viridiplantae; Streptophyta; Embryophyta; Tracheophyta;
Spermatophyta; Magnoliophyta; eudicotyledons; core eudicotyledons;
Caryophyllales; Cactaceae; Opuntioideae; Austrocylindropuntia.
REFERENCE 1 (bases 1 to 892)
AUTHORS Butterworth,C.A. and …我知道这应该很简单,但我不能让它工作.我的文件看起来像这样
>c12345|random info goes here that I want to delete
AAAAATTTTTTTTCCCC
>c45678| more | random info| here
GGGGGGGGGGG
而我想要做的就是让它变得更加简单,所以它可能看起来像这样
>seq1 [organism=human]
AAAAATTTTTTTTCCCC
>seq2 [organism=human]
GGGGGGGGGGGG
>seq3 [organism=human]
etc....
我知道一旦我通过执行以下操作获得索引部分,我就可以轻松追加该常量:
sed '/^>/ s/$/\[organism-human]/g'
但是如何构建该索引呢？
我在R中有一个数据框,如下所示:
Genes snps X0 X1 X2 X3
2 WASH7P 1_14677 0 2 2 2
3 WASH7P 1_14684 0 1 2 0
4 WASH7P 1_14685 0 0 0 0
是否有可能进行条件替换,如果int 2的列X0-X3的频率> 0.5,则将0替换为2,将2替换为0？这样新的数据帧是:
Genes snps X0 X1 X2 X3
2 WASH7P 1_14677 2 0 0 0
3 WASH7P 1_14684 0 1 2 0
4 WASH7P 1_14685 0 0 0 0
提前致谢!
我想处理多个床文件以查找重叠区域。我将我的数据集读取为数据帧，以及如何有效地并行扫描两个数据集以检测重叠区域发生的位置。我的方法是每次将数据框对象的每个单元格的峰值区域作为查询时，将间隔树中另一个数据框的所有行的峰值区域作为查询，然后搜索重叠区域。我很困惑如何在 R 中实现这一点。请帮助处理生物信息学中的床格式文件。感谢有人指出我如何做到这一点......
这是我想要实现的简单示例：
[1] chr1 [10171, 10226] * | MACS_peak_1 7.12
[2] chr1 [32698, 33079] * | MACS_peak_2 13.92
[3] chr1 [34757, 34794] * | MACS_peak_3 6.08
[4] chr1 [37786, 37833] * | MACS_peak_4 2.44
[5] chr1 [38449, 38484] * | MACS_peak_5 3.61
[6] chr1 [38584, 38838] * | MACS_peak_6 4.12
..
..
[] chrX [155191467, 155191508] * | MACS_peak_77948 3.80
[] chrX [155192786, 155192821] * | MACS_peak_77949 3.71
[] chrX [155206352, 155206433] * | MACS_peak_77950 3.81
[] chrX …我有一个像这样的字符串:
GATATATGGCACAGCAGTTGGATCCTTGAATGTCC[A/C]AGGTATATGTTATAGAAGCCTCGCAATTGTGTGTT
我意识到我可以使用替代这种特殊情况.但是,[X/Y]字符串之间的位置不同,括号之间的内容长度不同.
我想提取一下[X/Y].
我不想先解压缩vcf.gz文件,因为它是一个大文件而我的笔记本电脑没有空间.我试过做:
gunzip -c file.vcf.gz > bgzip -c > file.vcf.bgz
但它没有用.思考？
我有一个脚本,它获取查询ID的列表文件,并从uniprot中提取有机体和序列,代码运行良好,但它非常慢.我希望通过它处理大约400万个序列,但解析100个序列大约需要5分钟:
real 5m32.452s
user 0m0.651s
sys 0m0.135s
代码使用python的检索模块.我在网上看到我可以使用.session()属性,但是当我尝试这个时,我收到以下错误:
Traceback (most recent call last):
File "retrieve.py", line 14, in <module>
result = session.get(baseURL, payload)
TypeError: get() takes exactly 2 arguments (3 given)
代码列在这里:
import requests
baseURL = 'http://www.uniprot.org/uniprot/'
sample = open('sample.txt','r')
out = open('out','w')
for line in sample:
query = line.strip()
payload = {
'query': query,
'format':'tab',
'columns': 'id, entry_name, organism, sequence'
}
result = requests.get(baseURL, payload)
if result.ok:
out.write(query + '\t' + result.text[41:] + '\n')
输入格式示例:
EDP09046
ONI31767
ENSFALT00000002630
EAS32469 …在执行量化（kallisto/salmon）后，我正在使用一些快速的 R 脚本来 cbind 文件。问题是，我收到一个 R 错误，说我的输入文件长度不同，所以 cbind() 不起作用。
这显然不是这种情况，它们都是 16887 行（用 bash wc 检查），并且在没有snakenmake 的 R 中工作得很好。
还值得一提的是，我确实得到了随机数样本（~ 1 到 4）的输出。
这是R代码：
sample <- read.table(snakemake@input[[1]], sep = "\t", header = TRUE)
if (!file.exists(snakemake@params[[1]])) {
merge <- data.frame(sample)
} else {
merge1 <- read.table(snakemake@params[[1]], sep = "\t", header = TRUE)
merge <- cbind(merge1, sample[,2, drop=F])
}
write.table(merge, snakemake@params[[1]], sep = "\t", quote = F, row.names = F)
file.create(snakemake@output[[1]])
而蛇形规则：
rule merge_quantif:
input:
QUANTIF+"/{sample}/quantif.txt"
output:
QUANTIF+"/{sample}/merge.done"
params:
QUANTIF+"/all_sample_quantified.txt"
script:
"Tools/merge_quantif.R"
我的文件是这样的： …