标签: bioinformatics

这段代码我需要一些帮助.我知道应该递归的部分,或者至少我认为我这样做但不确定如何实现它.我正在尝试从对齐矩阵实现路径查找程序,该矩阵将找到多个路由回零值.例如,如果您执行我的代码并插入CGCA作为第一个序列,并将CACGTAT作为第二个序列插入,并将1,0和-1作为匹配,不匹配和差距分数.该程序给出了HDHHDD和对象的路径

CACGTAT

CGC - A-.

然而,除此之外,还有更多可能的路径和对象,除了我不知道有多少.我想要做的是让我的代码循环回到自身并找到其他路径和对齐,使用与第一次相同的代码,直到它用完可能的对齐.我在网上找到的最好的方法是递归,除了没有人能解释如何做到这一点.在这种情况下,应该有两个以上的路径和对话框HDDDHHD和CACGTAT,以及C - GCA-和.HDDDDHH,CACGTAT和--CGCA-.我只是不知道如何编写代码来执行此任务.

# Implementation of Needleman and Wunsch Algorithm

my($seq1, $len1, $seq2, $len2, $data, @matrix, $i, $j, $x, $y, $val1, $val2);
my($val3, $pathrow, $pathcol, $seq1loc, $seq2loc, $gapscore, $matchscore, $mismatchscore);

#first obtain the data from the user. 
print "Please enter the first sequence for comaprsion\n";
$seq1=<STDIN>;
chomp $seq1;

print "Please enter the second sequence for comparsion\n";
$seq2=<STDIN>;
chomp $seq2;


# adding extra characters so sequences align with matrix
# saves some calculations later on
$seq1 …

我已经创建了一个程序,可以读取可以产生互补链的DNA序列,进一步转化为mRNA.但是,我必须为该DNA找到最长的开放阅读框架.我编写了一些东西,但是当它打印出来时,我得不到答案.救命？

这就是我的意思;

# Search for the longest open reading frame for this DNA.
print "\nHere is the largest ORF, from 5' to 3':\n" ;
local $_ = $RNA_seq ;
while ( / AUG /g ) {
    my $start = pos () - 2 ;
    if ( / UGA|UAA|UAG /g ) {
        my $stop = pos ;
        $gene = substr ( $_ , $start - 1 , $stop - $start + 1 ), $/ ;
        print "$gene" ;
    }
}

# …

我有以下fasta文件:

'>gi|277456704|dbj|ID_P|Gene name LLL
MDGFAGSLDDSISAASTSDVQDRLSALESRVQQQEDEITVLKAALADVLRRLAISEDHVASVKKSVSSKV
YRRKHQELQAMQMELQSPEYKLSKLRTSTIMTDYNPNYCFAGKTSSISDLKEVPRKNITLIRGLGHGAFG
EVYEGQVSGMPNDPSPLQVAVKTLPEVCSEQDELDFLMEALIISKFNHQNIVRCIGVSLQSLPRFILLEL
MAGGDLKSFLRETRPRPSQPSSLAMLDLLHVARDIACGCQYLEENHFIHRDIAARNCLLTCPGPGRVAKI
GDFGMARDIYRASYYRKGGCAMLPVKWMPPEAFMEGIFTSKTDTWSFGVLLWEIFSLGYMPYPSKSNQEV
LEFVTSGGRMDPPKNCPGPVYRIMTQCWQHQPEDRPNFAIILERIEYCTQDPDVINTALPIEYGPLVEEE

'>gi|27704|dbj|ID_Y|Gene name JJJ
MDGFAGSLDDSISAASTSDVQDRLSALESRVQQQEDEITVLKAALADVLRRLAISEDHVASVKKSVSSKG
SELRGGYGDPGRLPVGSGLCSASRARLPGHVAADHPPAVYRRKHQELQAMQMELQSPEYKLSKLRTSTIM
TDYNPNYCFAGKTSSISDLKEVPRKNITLIRGLGHGAFGEVYEGQVSGMPNDPSPLQVAVKTLPEVCSEQ
DELDFLMEALIISKFNHQNIVRCIGVSLQSLPRFILLELMAGGDLKSFLRETRPRPSQPSSLAMLDLLHV
ARDIACGCQYLEENHFIHRDIAARNCLLTCPGPGRVAKIGDFGMARDIYRASYYRKGGCAMLPVKWMPPE

'>gi|2097704|dbj|ID_X|Gene name X
MDGFAGSLDDSISAASTSDVQDRLSALESRVQQQEDEITVLKAALADVLRRLAISEDHVASVKKSVSSKG
QPSPRAVIPMSCITNGSGANRKPSHTSAVSIAGKETLSSAAKSGTEKKKEKPQGQREKKEESHSNDQSPQ
IRASPSPQPSSQPLQIHRQTPESKNATPTKSIKRPSPAEKSHNSWENSDDSRNKLSKIPSTPKLIPKVTK
TADKHKDVIINQEGEYIKMFMRGRPITMFIPSDVDNYDDIRTELPPEKLKLEWAYGYRGKDCRANVYLLP
TGEIVYFIASVVVLFNYEERTQRHYLGHTDCVKCLAIHPDKIRIATGQIAGVDKDGRPLQPHVRVWDSVT
LSTLQIIGLGTFERGVGCLDFSKADSGVHLCVIDDSNEHMLTVWDWQRKAKGAEIKTTNEVVLAVEFHPT

我想循环通过FASTA,将蛋白质序列分解为它遇到的所有'R',这将生成肽,然后将肽进行爆破.从blastp获取结果并将blastp结果存储在fasta文件中每个蛋白质ID的单独文件中.我并不特别关注使用何种语言.我想学习如何做到这一点,以便我可以在它上面构建更多功能.谢谢!

我想在命令提示符下创建一个命令,通过输入'fetch ID.pdb'轻松下载蛋白质数据库文件.这些文件是从www.pdb.org/pdb/files/ID.pdb中检索的.

我尝试添加.bash_profiles;

    alias fetch ='wget www.pdb.org/pdb/files/'

我以为我可以输入'fetch ID.pdb'它会起作用,但似乎别名不能那样工作.有什么建议？如果我只能通过输入"获取ID"来检索蛋白质,那就更好了!

谢谢.

以下文件是配对端fastq文件的两个配对,我想根据它们的长度分隔每个fastq.

mate1.fq:

@SRR127.1
TGGTTATGATGTTTGTGTAGGAATAGAAATTTTGATTAAGATATTAGTGAAATTTGAATGTAGTTTATTTGGAAGTTATGGAGAGTTTATATTGTATTTATGTTTATTGTTGTAGATTTATATTTATGTGTATATATTAGTTTTTTTGTGT
+
ABAAAF4FFFFFGGGGGGFFGGFGHGFGHHHHHGGCFFGHHHHH5FDBED55DGGFEGFHHHGBHDDHHHFF3AB3FFG5CBGBEF5BD5DGFEGHFAGAFEDGHGFHHGHGEFFGFGGHFEGHHFHGBEBGHHHHGHBHHFHHGGFGHH2
@SRR127.2
TATGGTAAGAAAATTGAAAATTATAAAAAATGAAAAATGTTTATTTGATGATTTGAAAAATGATGAAATTATTGAAAAATGTGAAAAATGAGAAATGTATATTGTAGGATTTGGAATATGGTGAGATAAATGAAAATTATAGTAAATG
+
AABAA5@D4@5CFFCA55FFGGHDGFHFFCC45DGFA2FA5DD55AAAA55DDBDEDDBGGFF5BA5DDABF5D5B5FF1ADFB5EDGHFG5@BFBD55D5FFB@@5@GBGEFBGHHGB@DBBFHFBDG3B43FFH@FGFHH?FHHHH

mate2.fq:

@SRR127.1
ACCTATAAAAAAACCATATCAATAACTATAAAATCTTTATAAAATCCCACCCAATTAAAAAAAAATAAATTAATACATATAAAACCTTAAACACATAAAACATAATCACATACTATATAAACAATTACTATCACTACTAAACACCTAATA
+
>AA?AF13B@D@1EFCGGGFFG3EBGHHHBB2FGHHGHGFDGHHDFEGFHGGGHG1FFF1GGCGGGBGHHHHHFHHHHFHEGGFHF0BD1FGHHAGEGHFHHHFGGFHHGHHHFHHGGFHBGHFED1FBGFGFHDGHGHFGG1GB0GFHH
@SRR127.2
CTATTTCTCATTTTTTTATAATTTTCAATTCTCTTACCATATTCCACATCCTACACTAAACATTTCTAAATTTTCCACCTTTTTCTATTTTTCTCACCATATTTCATATCCTAAAAAACATATTCCTCATTTACTATAATTTTCAATTATC
+
11>>AFFDFF3@FFF?EFFGFBGHFDFA33D2FF2GGHFE12DD221AF1F1E1BG1GGBFBGGEGHDAABGAGDFABGG1BBDF12A2@2BG@2@DEFFF2B2@2222BB2211FGEE/11@22B2>1B22F2>GBGBD22BGD2>2B22

我编写了以下代码来执行此操作,但我只得到第二个文件(mate2.fq)的奇怪错误,而他们两个也有151 bp读取.

#!/usr/bin/perl

use strict;
use warnings;

my @fh;

my $file_name = $ARGV[0];
my $infile    = $ARGV[1];

#convert every 4-line fastq to 1-line
open(FH, "cat '$infile' | awk '{printf \"%s%s\",\$0,(NR%4?FS:RS)}' | ");

while (<FH>) {
  chomp;

  my @line = split(/\s+/, $_);
  my $len  = length($line[1]);

  if ($len >= 100) {

    #print $len,"\n",$_,"\n";
    push @fh, $len;

    if (not …

所以我试图转换一组矩阵格式的数据,就像这样

  V1  V2  V3  V4
1 12  32  12  54
2 43  34  23  24
3 43  23  65  54

其中1-3是seqID,V1-4是ID

我希望以一种更容易让SPSS理解的格式,如下所示:

ID  seqID  Value
1      1   12
2      1   32
3      1   12
4      1   54
1      2   43

等等

R中最好的方法是什么？它应该是一些简单的文本操作,我尝试制作一个循环来解决它(例如通过提取列,然后再次绑定值)但不知何故循环语法仍然逃避我!帮助赞赏!

我需要perl中的正则表达式,它在开始时与ATG匹配,并以TAG,TAA或TGA结束.这是我到目前为止的代码.

my $sequence = 'AATGGTTTCTCCCATCTCTCCATCGGCATAAAAATACAGAATGATCTAACGAA';

while($sequence =~ ____) {
 print $1;
}

我认为这对你来说是一个简单的问题,因为我是python3的初学者.打印fasta文件的标题时,它包含括号.我怎么能删除它们？

import sys
from Bio import Entrez
from Bio import SeqIO


#define email for entrez login
db           = "nuccore"
Entrez.email = "someone@email.com"

#load accessions from arguments
if len(sys.argv[1:]) > 1:
  accs = sys.argv[1:]
else: #load accesions from stdin  
  accs = [ l.strip() for l in sys.stdin if l.strip() ]
#fetch
sys.stderr.write( "Fetching %s entries from GenBank: %s\n" % (len(accs), ", ".join(accs[:10])))
for i,acc in enumerate(accs):
  try:
    sys.stderr.write( " %9i %s          \r" % (i+1,acc))  
    handle = Entrez.efetch(db=db, rettype="fasta", id=acc)
    seq_record = …

这个问题应该很简单,但我不知道如何在R中完成.基本上我有一个双列文件,其中,第一列表示脚手架ID,第二列是SNP的位置那个特别的脚手架.

id    POS
0001  38
0001  46
0001  50
0002  17
0002  23
0002  46
0003  13
0003  19
0003  38
...   ...

我想创建一个新列,其中第一个支架的SNP位置(NEW_POS)是相同的,但对于第二个支架(以及之后),NEW_POS将是添加前一个支架的最后一个POS值的结果第二个支架中每个POS的支架(即50 + 17,50 + 23,50 + 46,......).对于第三个脚手架(96 + 13,96 + 19,96 + 38,......)等.正如你在这里看到的:

id    POS  NEW_POS
0001  38   38
0001  46   46
0001  50   50
0002  17   67
0002  23   73
0002  46   96
0003  13   109
0003  19   115
0003  38   134
...   ...  ...

如何可视化两个序列的完全比对？

library(Biostrings)
s1 <-DNAString("ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAGAAGACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGTTTTCAAG")
s2 <-DNAString("GTTTCACTACTTCCTTTCGGGTAAGTAAATATATGTTTCACTACTTCCTTTCGGGTAAGTGTTTCACTACTTCCTTTCGGGTAAGTAAATATATAAATATATAAAAATATAATTTTCATCAAATATATAAATATATAAAAATATAATTTTCATCAAATATATAAAAATATAATTTTCATC")
pairwiseAlignment(s1,s2)

输出：

Global PairwiseAlignmentsSingleSubject (1 of 1)
pattern: [1] ACTTCACCAGCTCCCTGGCGGTAAGTTGATCAAAGGAAACGCAAAGT--TTTCAC---...CTTCACCAGCTCCCTGGCGGTAAGTTG-ATCAAAGG---AAACGCAAAGTTTTCAAG 
subject: [1] GTTTCACTACTTCCTTTCGGGTAAGTAAAT-ATATGTTTCACTACTTCCTTTCGGGTA...TATATAAATATATAAAAATATAATTTTCATCAAATATATAAAAATATAATTTTCATC 
score: -394.7115 

在这里，仅显示了部分对齐？您知道绘制或打印路线的任何现有功能吗？