连接外显子序列并在其间插入 Ns
Concatenate exon sequences and inserting Ns in between
所以我想用带有外显子序列的 FASTA 文件执行某种复杂的任务,例如:
>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"、外显子编号和字符串长度组成。
我想要实现的是:
>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。比如说,在第一种情况下,引入一串31个Ns在外显子 1 和 2 之间。而在第二种情况下,在外显子 1 和 2 之间引入一串 31 Ns,在外显子 2 和 3 之间引入一串 16 Ns。
这基本上是我的主要任务,也是棘手的任务。有谁知道或提出基于 Python、BASH、AWK、R、Perl 或类似的解决方案?
我一直在尝试解决这个问题,但我无法找到一个好的通用解决方案来在合并它们时循环相同的转录本...
非常感谢。
Perl 解决方案:
#!/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) = (, );
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 及其部分(以及前一个外显子)存储在变量中。如果遇到 header,我们检查是否有新的转录本开始。如果是,我们可以忘记前面的外显子,否则我们打印 header 包含转录本,前面的部分和实际的部分。
如果我们读取外显子行,如果有一个外显子要打印(这意味着我们不在转录本的第一部分),我们会检查新旧外显子的长度,select越长的,输出外显子和N's。然后,我们记住外显子以进行进一步处理,如果同一转录本的另一部分跟随。
所以我想用带有外显子序列的 FASTA 文件执行某种复杂的任务,例如:
>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"、外显子编号和字符串长度组成。
我想要实现的是:
>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。比如说,在第一种情况下,引入一串31个Ns在外显子 1 和 2 之间。而在第二种情况下,在外显子 1 和 2 之间引入一串 31 Ns,在外显子 2 和 3 之间引入一串 16 Ns。
这基本上是我的主要任务,也是棘手的任务。有谁知道或提出基于 Python、BASH、AWK、R、Perl 或类似的解决方案?
我一直在尝试解决这个问题,但我无法找到一个好的通用解决方案来在合并它们时循环相同的转录本...
非常感谢。
Perl 解决方案:
#!/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) = (, );
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 及其部分(以及前一个外显子)存储在变量中。如果遇到 header,我们检查是否有新的转录本开始。如果是,我们可以忘记前面的外显子,否则我们打印 header 包含转录本,前面的部分和实际的部分。
如果我们读取外显子行,如果有一个外显子要打印(这意味着我们不在转录本的第一部分),我们会检查新旧外显子的长度,select越长的,输出外显子和N's。然后,我们记住外显子以进行进一步处理,如果同一转录本的另一部分跟随。