首页
学习
活动
专区
圈层
工具
发布
首页
学习
活动
专区
圈层
工具
MCP广场
社区首页 >问答首页 >连接外显子序列并在其间插入Ns

连接外显子序列并在其间插入Ns
EN

Stack Overflow用户
提问于 2019-07-04 01:36:35
回答 1查看 81关注 0票数 1

所以我想用一个带有外显子序列的FASTA文件来执行一个复杂的任务,例如:

代码语言:javascript
运行
复制
>MSTRG.1.1_1_30
AAAACGGAGGACCAAGCCGTTTGCCGTACG
>MSTRG.1.1_2_20
CCCGAAGGGCGTTAGTGAGC
>MSTRG.2.1_1_30
ACGGGAGCGTTGTCGACCGGTTGCGAGCGT
>MSTRG.2.1_2_10
AACGGGACCT
>MSTRG.2.1_3_15
AACGTTTGCGTCTCT

可以看到,我有两个文本,称为MSTRG.1.1和MSTRG.2.1。第一个文本有两个外显子,第一个被称为MSTRG.1.1_1_30,长度为30个字母。第二个外显子(片段)有20个字母。相反,另一个转录本有三个外显子,每个外显子有30,10和15个字母。

有更多的转录本和超过3个外显子,最多20个或多或少。

正如您所看到的,字母序列的标识符由抄本的名称"MSTRG.X.X“、外显子的编号和字符串的长度组成。

我想实现的目标如下:

代码语言:javascript
运行
复制
>MSTRG.1.1_1_2
AAAACGGAGGACCAAGCCGTTTGCCGTACGNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNCCCGAAGGGCGTTAGTGAGC
>MSTRG.2.1_1_2
ACGGGAGCGTTGTCGACCGGTTGCGAGCGTNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNAACGGGACCT
>MSTRG.2.1_2_3
AACGGGACCTNNNNNNNNNNNNNNNNAACGTTTGCGTCTCT

一些解释:

我想把所有连续的外显子组合在同一个转录本中,所以对于转录本1,因为有两个外显子,所以只有一个合并的序列是由外显子1和2组合产生的。对于转录本2,如果存在更多的外显子,结果将是两个组合,外显子1和2,以及外显子2和3,依此类推,从而得到n-1个组合,其中n是外显子的数量。

除此之外,我想在两个组合外显子之间引入一串Ns,长度等于合并的外显子对的最长外显子+ 1。例如,在第一种情况下,在外显子1和2之间引入一串31 Ns。而在第二种情况下,在外显子1和2之间引入一串31 Ns,在外显子2和3之间引入一串16 Ns。

这基本上是我的主要任务,也是棘手的任务。有没有人知道或提出了基于Python、BASH、AWK、R、Perl或类似的解决方案?

我一直在尝试解决这个问题,但我无法找到一个很好的通用解决方案来循环相同的脚本,同时合并它们……

非常感谢。

EN

回答 1

Stack Overflow用户

回答已采纳

发布于 2019-07-04 02:00:02

Perl解决方案:

代码语言:javascript
运行
复制
#!/usr/bin/perl
use warnings;
use strict;
use feature qw{ say };

my $transcript;
my $part;
my $exon;
while (<>) {
    chomp;
    if (/^>(.*?\.[0-9]+\.[0-9]+)_([0-9]+)/) {
        my ($new_transcript, $new_part) = ($1, $2);
        if ($transcript && $transcript eq $new_transcript) {
            say ">$transcript\_$part\_$new_part";

        } else {
            undef $exon;
        }

        $transcript = $new_transcript;
        $part = $new_part;

    } else {
        my $new_exon = $_;
        if ($exon) {
            my $max_length = length $exon;
            $max_length = length $new_exon if length $new_exon > $max_length;
            say $exon, 'N' x ++$max_length, $new_exon;
        }
        $exon = $new_exon;
    }
}

它逐行读取文件。最后处理的转录id及其部分(以及前一个外显子)存储在变量中。如果我们遇到一个标题,我们检查是否启动了一个新的脚本。如果是,我们可以忘记前一个外显子,否则我们打印包含抄本、前一部分和实际部分的标题。

如果我们阅读外显子行,如果有一个外显子需要打印(这意味着我们不是在转录的第一部分),我们检查旧的和新的外显子的长度,选择较长的外显子,并输出外显子和N's。然后,我们记住外显子以便进一步处理,如果后面跟着相同转录的另一部分。

票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/56875399

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档