我正试图在http://rosalind.info/problems/iprb/解决问题
鉴于:三个正整数
k,m和n代表含有k+m+n生物的群体:k个体是一个因子的纯合优势,m是杂合的,并且n是纯合的隐性.返回:两个随机选择的交配生物将产生具有显性等位基因的个体(从而显示显性表型)的概率.假设任何两种生物都可以交配.
我的解决方案适用于样本,但不适用于生成的任何问题.在进一步研究之后,似乎我应该找到随机选择任何一个生物体的概率,找到选择第二个生物体的概率,然后是产生后代与显性等位基因的配对概率.
我的问题是:我的代码在下面找到了什么概率？是否找到了所有可能交配的具有显性等位基因的后代的百分比 - 因此,如果对所有对进行测试,我的代码是否正在解决具有显性等位基因的后代的百分比？
f = open('rosalind_iprb.txt', 'r')
r = f.read()
s = r.split()
############# k = # homozygotes dominant, m = #heterozygotes, n = # homozygotes recessive
k = float(s[0])
m = float(s[1])
n = float(s[2])
############# Counts for pairing between each group and within groups
k_k = 0
k_m = 0
k_n = 0
m_m = 0 …我有不同组织的数据,如此
tissueA tissueB tissueC
gene1 4.5 6.2 5.8
gene2 3.2 4.7 6.6
我想计算一个汇总统计量
x = ? [1-log2(i,j)/log2(i,max)]/n-1
其中n是组织的数量(这里是3),(i,max)是n个组织中基因i的最高值(即对于gene1,它是6.2).
因为我必须为每个基因的每个组织做这个计算(因为总和从j到n,并且j = 1),然后得到它的总和
我写了一个for循环
for (i in seq_along(x) {
my.max <- max(x[,i])
my.statistic <- (1-log2(x[,i]/log2[my.max])
my.sum <- sum(my.statistic)
my.answer <- my.sum/2 #(n-1=3-1=2)
但是我不知道如何为每一行应用这个for循环,通常我会编写一个函数并且只执行(apply,1,function(x))但是我不确定如何将for循环转换为函数.
例如,对于gene1的预期输出,它就是
(1-log2(4.5)/log2(6.2))/2 + (1-log2(5.8)/log2(6.2))/2 =0.1060983
我想从multifasta文件中提取与单独的ID列表给出的ID匹配的序列.
FASTA文件seq.fasta:
>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11605
TTCAGCAAGCCGAGTCCTGCGTCGAGAGTTCAAGTC
CCTGTTCGGGCGCCACTGCTAG
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
>7P58X:01334:11635
TTCAGCAAGCCGAGTCCTGCGTCGAGAGATCGCTTT
CAAGTCCCTGTTCGGGCGCCACTGCGGGTCTGTGTC
GAGCG
>7P58X:01336:11621
ACGCTCGACACAGACCTTTAGTCAGTGTGGAAATCT
CTAGCAGTAGAGGAGATCTCCTCGACGCAGGACT
ID文件id.txt:
7P58X:01332:11636
7P58X:01334:11613
我想获取fasta文件,只包含与id.txt文件中的ID匹配的序列:
>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
我真的很喜欢我在这里和这里的答案中找到的awk方法,但是那里给出的代码仍然不能完美地用于我给出的例子.原因如下:
(1)
awk -v seq="7P58X:01332:11636" -v RS='>' '$1 == seq {print RS $0}' seq.fasta
此代码适用于多行序列,但ID必须单独插入代码.
(2)
awk 'NR==FNR{n[">"$0];next} f{print f ORS $0;f=""} $0 in n{f=$0}' id.txt seq.fasta
此代码可以从id.txt文件中获取ID,但只返回多行序列的第一行.
我想好的方法是修改代码中的RS变量(2),但到目前为止我的所有尝试都失败了.请问,有人帮我吗？
在一篇关于生命科学标识符的论文中(参见LSID Tester,一种测试生命科学标识符解析服务的工具),Roderic DM Page博士写道:
鉴于LSID urn:lsid**:ubio.org**:namebank:11815,查询DNS以获取_lsid._tcp的SRV记录.ubio.org返回animalia.ubio.org:80作为ubio.org LSID服务的位置.
我了解到我可以使用unix上的host命令将_lsid._tcp.ubio.org链接到animalia.ubio.org:80 :
host -t srv _lsid._tcp.ubio.org
_lsid._tcp.ubio.org has SRV record 1 0 80 ANIMALIA.ubio.org
如何使用Java J2SE API执行此"DNS"操作(没有任何外部Java库,我想要一个轻量级的解决方案)？
谢谢
我被教授HMM并且给了这个家庭作业问题.我理解了它的一部分,但我不确定它是否正确.问题是:
考虑一个不同的游戏,经销商不会翻转硬币,而是使用标签1,2和3滚动三面模具.(尽量不要考虑三面模具的外观.)经销商已经两个装载的骰子D1和D2.对于每个骰子Di,滚动数字i的概率是1/2,并且其他两个结果中的每一个的概率是1/4.在每个回合,经销商必须决定是(1)保持相同的骰子,(2)切换到另一个骰子,或(3)结束游戏.他选择(1)概率为1/2,其他每个概率为1/4.开始时,经销商以相同的概率选择两个骰子中的一个.
为这种情况提供HMM.指定字母表,状态,转换概率和发射概率.包括开始状态开始,并假设HMM以状态开始以概率1开始.还包括结束状态结束.
假设您观察到以下的模具辊序列:1 1 2 1 2 2.找到最能说明辊子顺序的状态序列.这个序列的概率是多少？通过完成Viterbi表找到答案.在单元格中包含回溯箭头,以便您可以追溯状态序列.以下某些事实可能有用:
log2(0)=-∞log2
(1/4)= -2
log2(1/2)= -1
log2(1)= 0- 对于这种模具辊序列,实际上存在两种最佳状态序列.另一个国家的序列是什么？
如果我对第一部分没有错,我必须做类似这里的事情http://en.wikipedia.org/wiki/Hidden_Markov_model#A_concrete_example但是我没有得到假设以概率1开始的东西.
另外,我不知道在问题的第二部分我要为维特比表做些什么.如果任何身体可以给我一些提示或线索,那将是伟大的.
我正在尝试Bio::DB::Sam在远程服务器上的主目录上安装perl模块.
我下载了模块,解压缩了文件,然后运行:
perl Build.pl prefix=~/local
接下来会发生什么:
This module requires samtools 0.1.10 or higher (samtools.sourceforge.net).
Please enter the location of the bam.h and compiled libbam.a files: **/some_places/samtools-0.1.19**
Found /some_places/samtools-0.1.19/bam.h and /some_places/samtools-0.1.19/libbam.a.
Created MYMETA.yml and MYMETA.json
Creating new 'Build' script for 'Bio-SamTools' version '1.39'
接下来当我尝试运行时:
./Build
这就是我得到的:
Building Bio-SamTools
gcc -shared -O2 -g -pipe -Wall -Wp,-D_FORTIFY_SOURCE=2 -fexceptions -fstack-protector --param=ssp-buffer-size=4 -m64 -mtune=generic -o blib/arch/auto/Bio/DB/Sam/Sam.so lib/Bio/DB/Sam.o c_bin/bam2bedgraph.o -L/some_places/samtools-0.1.19 -lbam -lpthread -lz
/usr/bin/ld: /some_places/samtools-0.1.19/libbam.a(bgzf.o): relocation R_X86_64_32 against `.rodata.str1.1' can not be used when …科学名称通常由3个信息组成:属,种上皮和作者.一个简单的例子如下:
Acanthus ilicifolius L.
简单.然而,当我们必须处理杂交种,亚种/品种/形态,几个作者和其他不一致时,问题变得更加复杂.在这些情况下,物种名称可能如下所示:
比照 穿心莲(Burm.f.)墙.前Nees
或这个:
Ipomoea pes-caprae(L.)DC.亚种.brasiliensis(L.)Ooststr.f
我正试图找到一种可靠的解构这些名字的方法.如果if/else语句我可以使用吨写一些hackish代码,但我正在寻找更优雅(和健壮)的东西.我想的是某种解析器的名称类似于解析数学表达式的计算器.不幸的是,我不是最复杂的程序员,我之前也没有写过真正的解析器,也不知道在这种情况下它是否有意义,因为科学名称有很多变化.你认为解决这个问题的最佳方法是什么？首选语言是R,如果它更适合任务,也许也是Julia.
下午好,我试图计数的时间ACTG发生在DNA序列使用perl6.i字母都试过其他方式我只是试图让它以另一种方式完成的数量.以下是我提出的一些代码
use v6;
my $default-input = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC";
sub MAIN(Str $input = $default-input)
{
say "{bag($input.comb)<A C G T>}";
}
use v6;
my $default-input = "AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC";
sub MAIN($input = $default-input)
{
"{<A C G T>.map({ +$input.comb(/$_/) })}".say;
数据集样本
AGCTTTTCATTCTGACTGCAACGGGCAATATGTCTCTGTGTGGATTAAAAAAAGAGTGTCTGATAGCAGC
对于初学者我不熟悉生物信息学,特别是编程,但我已经构建了一个脚本,它将通过一个所谓的VCF文件(只包括个体,一个clumn =一个人),并使用搜索字符串来查找对于每个变体(系),个体是纯合的还是杂合的.
这个脚本起作用,至少在小子集上,但我知道它将所有东西都存储在内存中.我想在非常大的压缩文件(甚至整个基因组)上做这个,但我不知道如何将这个脚本转换成一个逐行完成所有操作的脚本(因为我想要计算整列,我只是不要看看如何解决这个问题.
因此,每个人的输出是5件事(总变体,数字纯合子,数字杂合子,以及同源和杂合子的比例).请参阅以下代码:
#!usr/bin/env python
import re
import gzip
subset_cols = 'subset_cols_chr18.vcf.gz'
#nuc_div = 'nuc_div_chr18.txt'
gz_infile = gzip.GzipFile(subset_cols, "r")
#gz_outfile = gzip.GzipFile(nuc_div, "w")
# make a dictionary of the header line for easy retrieval of elements later on
headers = gz_infile.readline().rstrip().split('\t')
print headers
column_dict = {}
for header in headers:
column_dict[header] = []
for line in gz_infile:
columns = line.rstrip().split('\t')
for i in range(len(columns)):
c_header=headers[i]
column_dict[c_header].append(columns[i])
#print column_dict
for key in column_dict:
number_homozygotes = 0
number_heterozygotes = 0
for …我试图.bed通过识别前两列chr并start遵循此来合并多个文件，
但是，我想知道如何使文件名成为新添加的列名。
$cat combineFWPS_02.sh
BEGIN {
for (k=1; k<ARGC; ++k)
s = s " " 0
}
FNR == 1 {
++ARGIND
}
{
key=$1 OFS $2
if (!(key in map))
map[key] = s
split(map[key], a)
a[ARGIND] = $3
v = ""
for (k=1; k<ARGC; ++k)
v = v " " a[k]
map[key]=v
}
END {
for (k in map)
print k map[k]
}
$cat comRwps_02.sh
awkCOM="~/scripts/combineFWPS_02.sh"
## Run the jobs
time awk …