Chinaunix

标题: 20积分悬赏，求个perl程序，新手，写了好几个都无果。 [打印本页]

作者: yanglc2013 时间: 2013-03-27 19:13
标题: 20积分悬赏，求个perl程序，新手，写了好几个都无果。
本帖最后由 yanglc2013 于 2013-04-03 15:52 编辑

其中ref图中的数据为库，在同一条染色体上（两个图片中蓝色区域必须先相同），如果lnc图中红色矩形内的区域包含在ref图中红色区域内输出文件1（contain），

两个红色矩形框内的数据有交集输出文件2（overlap），无交集输出文件3（gapped）.如果lnc区间比ref大，就是ref图中红色区域包含在lnc图中红色区域输出文件4（reserve）.

希望输出的格式：（比如下面是overlap，不仅要输出lnc的数据，还有在同一行加上ref对应的数据。）

lnc-C1orf227-1:2 chr1 213028819 213031452 - NM_014053 chr1 213031596 213072705 + 2210

后面的高手只是输出了lnc的数据。

数据文件.rar

920.59 KB, 下载次数: 38

作者: picbhan 时间: 2013-03-27 19:13

#!/usr/bin/perl
use strict;
use warnings;
my (%ref, %data);
open my $in , '<', 'refGene.txt'
or die "Can't open 'refGene.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$ref{$t[1]}}, [@t];
}
close $in;
open $in , '<', 'lncipedia_location.txt'
or die "Can't open 'lncipedia_location.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$data{$t[1]}}, [@t];
}
close $in;
open my $include , '>', 'include.txt'
or die "Can't create 'include.txt' for writing. $!\n";
open my $non_include , '>', 'non_include.txt'
or die "Can't create 'non_include.txt' for writing. $!\n";
open my $overlap , '>', 'overlap.txt'
or die "Can't create 'overlap.txt' for writing. $!\n";
while (my ($chr, $v) = each %ref) {
my @ref_data = sort {$a->[2] <=> $b->[2]} @$v;
next unless exists $data{$chr};
my @inc_data = sort {$a->[2] <=> $b->[2]} @{$data{$chr}};
for (@ref_data) {
my @rd = @$_;
last unless @inc_data;
while (1) {
last unless @inc_data;
my @id = @{+shift @inc_data};
if ($id[3] < $rd[2]) {
print {$non_include} join("\t", @id), "\n";
next;
}
if ($id[2] > $rd[3]) {
unshift @inc_data, [@id];
last;
}
if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
print {$include} join("\t", @id), "\n";
}
else {
print {$overlap} join("\t", @id), "\n";
}
}
}
delete $data{$chr};
}
while (my ($chr, $v) = each %data) {
my @inc_data = @$v;
for (@inc_data) {
print {$non_include} join("\t", @$_), "\n";
}
}
close $include;
close $non_include;
close $overlap;

复制代码

作者: Perlvim 时间: 2013-03-27 20:44
本帖最后由 Perlvim 于 2013-03-27 20:48 编辑

是测试DNA相似度吗？

1. refGene 文件中的
NR_026818 chr1 34610 36081 - 585
NR_033908 chr1 700244 714068 - 590
6个字段，是否只使用第2个到第4个，其他没有用？

2. refGene 文件中第二个字段相同的，是否要合并范围？计算并集？

3. lnc 中的数据，计算在哪个范围，是否有交集，是在 ref 那个对应字段全部范围内计算吗？

作者: mcshell 时间: 2013-03-27 20:50
不懂生物真苦逼{:3_186:}

作者: yanglc2013 时间: 2013-03-27 22:02
回复 2# Perlvim
一个染色体可以分好几段，虽然你说的有重复，但是其起始位置和终止位置不一样，所以不要合并啊。这是求gene的location信息。判断lnc中的染色体（chr）其起始位置和终止位置是不是包含在ref中的起始位置和终止位置中。就是你指的第三列（起始位置）和第四列（终止位置）。

作者: Perlvim 时间: 2013-03-27 22:41
回复 5# yanglc2013

意思是一个染色体中一个位置有很多片段，要每个片段都查找，是否有交集吗？

作者: yanglc2013 时间: 2013-03-27 23:12
回复 6# Perlvim
对啊。就是比对染色体，如果lnc中的那一段染色体在ref中就输出，不在输出待另一个文件，有一部分在又输出一个文件。

作者: dahe_1984 时间: 2013-03-27 23:20
数据好大，机器都快跑死了，你这个够牛的，不好做啊。

作者: blackfur 时间: 2013-03-28 00:53
回复 1# yanglc2013

神马叫‘落在’，是两个数值相同的意思吗，生物不大懂，建议楼主先扫盲

作者: dahe_1984 时间: 2013-03-28 08:32

blackfur 发表于 2013-03-28 00:53
回复 1# yanglc2013

落在就是包含呗。比如 4..6 落在 1..10中

作者: yanglc2013 时间: 2013-03-28 08:58
回复 9# blackfur

意思就是ref中的红色区间包含lnc那个红色区间啊。其实就是两个相同的染色体随机切割后。有的会没有交集，有的有，有的一条包含在另一条内。

作者: yanglc2013 时间: 2013-03-28 09:02
回复 9# blackfur
楼下正解啊

作者: picbhan 时间: 2013-03-28 09:27
昨晚写的那个程序忘了三行，结果漏掉了一点没有交集的行。修改如下：

#!/usr/bin/perl
use strict;
use warnings;
my (%ref, %data);
open my $in , '<', 'refGene.txt'
or die "Can't open 'refGene.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$ref{$t[1]}}, [@t];
}
close $in;
open $in , '<', 'lncipedia_location.txt'
or die "Can't open 'lncipedia_location.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$data{$t[1]}}, [@t];
}
close $in;
open my $include , '>', 'include.txt'
or die "Can't create 'include.txt' for writing. $!\n";
open my $non_include , '>', 'non_include.txt'
or die "Can't create 'non_include.txt' for writing. $!\n";
open my $overlap , '>', 'overlap.txt'
or die "Can't create 'overlap.txt' for writing. $!\n";
while (my ($chr, $v) = each %ref) {
my @ref_data = sort {$a->[2] <=> $b->[2]} @$v;
next unless exists $data{$chr};
my @inc_data = sort {$a->[2] <=> $b->[2]} @{$data{$chr}};
for (@ref_data) {
my @rd = @$_;
last unless @inc_data;
while (1) {
last unless @inc_data;
my @id = @{+shift @inc_data};
if ($id[3] < $rd[2]) {
print {$non_include} join("\t", @id), "\n";
next;
}
if ($id[2] > $rd[3]) {
unshift @inc_data, [@id];
last;
}
if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
print {$include} join("\t", @id), "\n";
}
else {
print {$overlap} join("\t", @id), "\n";
}
}
}
for (@inc_data) {
print {$non_include} join("\t", @$_), "\n";
}
delete $data{$chr};
}
while (my ($chr, $v) = each %data) {
my @inc_data = @$v;
for (@inc_data) {
print {$non_include} join("\t", @$_), "\n";
}
}
close $include;
close $non_include;
close $overlap;

复制代码

作者: picbhan 时间: 2013-03-28 09:28
回复 8# dahe_1984

为毛我跑的时候貌似一秒钟就完了，数据不大啊。

作者: yanglc2013 时间: 2013-03-28 09:33
回复 4# picbhan
这个程序好像无论ref还是lnc中重复的chr没有比对啊。那不能忽略，其起始和终止位置不一样，还是得单独比对。输出的文件中，可能三个文件中有相同的name，但是其区间不一样。

作者: picbhan 时间: 2013-03-28 09:40
回复 15# yanglc2013

我不明白你说的什么意思。我是按照你在帖子开头说的意思来写代码的，即只在inc和ref之间比较。所以你如果觉得不是这个意思，麻烦说清楚你到底要什么结果。从生物的角度说就行，你这每个区间是个基因？

作者: picbhan 时间: 2013-03-28 09:42
回复 15# yanglc2013

你说的name是指第一列吗？同一个文件里面相同染色体上的区间也要比较？

作者: climby 时间: 2013-03-28 10:14
如果让我迅速的解决这个问题，鉴于数据的工整性，我会将数据倒入数据库，利用数据库的查询功能，会非常迅速的得到想要的结果。
如果非要用Perl 来解决这个问题，我看了一下数据， ref文件有4万多条记录， lnc文件有3万多条记录，一次读进内存然后比较的话，对机器的要求还是比较高的。
我的处理方法是
1. 先按照第二列对文件进行分组处理。
读取ref 文件，按照第二列进行分组，生成一些小文件，放在ref目录下。
分组后的文件如下：

对lnc文件进行类似处理

2.然后对分组后的文件一对一进行对比分析。

这样处理起来，应该对机器的压力会小很多。由于我上午抽不出时间来完成所有代码，仅列出方法供大家参考

作者: rubyish 时间: 2013-03-28 11:18
本帖最后由 rubyish 于 2013-03-31 12:50 编辑

是否:

Delete

有错误

复制代码

作者: yanglc2013 时间: 2013-03-28 11:39
回复 18# climby
这个方法我也试过，但是比较麻烦，因为我后面还有很多这样的数据。如果每个都得细分，那很耗时间的。

作者: climby 时间: 2013-03-28 11:49
回复 20# yanglc2013

反正是用脚本划分文件，效率应该不会低，这个4万多条记录的文件，划分小文件，差不多1秒就完成了。这种划分小文件是为后面的对比做准备，提高后面对比的效率。

作者: yanglc2013 时间: 2013-03-28 12:00
回复 16# picbhan

区间是基因的起始和终止位置。就是比如，lnc中有很多chr1，ref中也有很多chr1，lnc中的每个chr1都要和ref的chr1分别比较。

QQ截图20130328115513.png (4.9 KB, 下载次数: 32)

作者: yanglc2013 时间: 2013-03-28 12:02
回复 21# climby
谢谢啊，下次我用这个方法。

作者: picbhan 时间: 2013-03-28 12:29
回复 22# yanglc2013

好吧，明白你说的了。但是显然你没看懂程序里的算法。因为你要的只是找出inc中每条染色体上每个区间是否与ref中相同染色体上的区间有交集（完全包含，部分交集，不包含），所以在算法上不需要在两个文件间一一比较，这样很费时间。我在读入文件的时候先按染色体分类，完了用每个区间的起始位置排序了，这样就不需要全部比较。就好像如果你要找出 2 在数组 1， 5， 7， 19， 21 中的位置，只需要找出大于和小于2的临界位置就行了，不需要和5以后的数据比较。如果还不理解我再给你仔细分析算法吧。

作者: yanglc2013 时间: 2013-03-28 13:23
回复 24# picbhan
谢了啊大哥。我测试了下，你的程序可以达到要求，但是还有一种情况就是 lnc中的范围超出了ref中的范围，就是，lnc的范围包含ref的范围，这种情况能帮我输出到另一个文件吗。您的结果中把ref包含lnc和lnc包含ref都输入到了同一个文件，我试着改了下，但没改好。

作者: yanglc2013 时间: 2013-03-28 13:25
回复 24# picbhan
您的思路确实很厉害啊，

作者: picbhan 时间: 2013-03-28 15:16
回复 25# yanglc2013

我想问一下就是你要的结果是不是只是inc这个文件的？或者是两个文件的？
比如inc的一个区间包含ref里的一个区间，我是输出ref的那个区间还是inc的？在原来的程序里我只输出inc里的区间，没输出任何在ref里的数据。在原来的程序里你要的最后一种情况的结果应该在overlap的数据里，而不是在include的文件里。其实原来的程序里有点bug。我写好再给你新的代码。

作者: yanglc2013 时间: 2013-03-28 16:13
回复 27# picbhan
还是输出lnc的。因为可能会出现这种情况。真谢谢你啊。看你的程序收货很大，排版好，思路也巧妙。

作者: yanglc2013 时间: 2013-03-28 16:14
回复 27# picbhan

ref只是一个数据库。

作者: yanglc2013 时间: 2013-03-28 16:52
本帖最后由 yanglc2013 于 2013-03-28 16:59 编辑

回复 2# picbhan
大哥，文件中不是还有 + -之分。表示正链和负链，同一个染色体（比如chr1），他可以是正链也可以是负链，如果我还得加上++，--，+-，-+（前一个正负号为lnc中，后一个为ref中）必须也相同这个条件（就是除了必须染色体（chr）相同，而且lnc的+和ref中的+，lnc的-和ref的-，lnc的+和ref的-，lnc的-和ref的+），输出的结果还是lnc中的数据，ref作为一个数据库。那应该在你的程序上怎么修改？

作者: picbhan 时间: 2013-03-28 18:25
回复 31# yanglc2013

按你的意思，你是要输出到16个文件中？四种情况：++，--，+-，-+，每种情况下四种条件：inc包含ref，ref包含inc，inc与ref部分overlap，没有overlap？
因为同一文件中相同染色体上不同基因之间本身就有overlap，这样可能导致inc中某一个peak可能出现两种甚至三种不同的条件，比如inc中某一个区间和ref中三个不同基因分别是inc包含ref，ref包含inc，inc与ref部分overlap的情况，那我到底该输出到哪个文件？你设的条件越多，inc中同一个区间出现不同情况的可能性越大。

作者: yanglc2013 时间: 2013-03-28 18:38
回复 19# rubyish
能加上注释吗

作者: yanglc2013 时间: 2013-03-28 18:40
回复 32# picbhan
不考虑正负啊。我只是想学学，应该怎么改程序。加上注释，我就可以看明白了。有些地方我不太懂。

作者: picbhan 时间: 2013-03-28 19:36
本帖最后由 picbhan 于 2013-03-28 19:58 编辑

#!/usr/bin/perl
use strict;
use warnings;
my (%ref, %data);
open my $in , '<', 'refGene.txt'
or die "Can't open 'refGene.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$ref{$t[1]}}, [@t];
}
close $in;
open $in , '<', 'lncipedia_location.txt'
or die "Can't open 'lncipedia_location.txt' for reading. $!\n";
while (<$in>) {
chomp;
my @t = split;
push @{$data{$t[1]}}, [@t];
}
close $in;
open my $ref_include_inc , '>', 'ref_include_inc.txt'
or die "Can't create 'ref_include_inc.txt' for writing. $!\n";
open my $inc_include_ref , '>', 'inc_include_ref.txt'
or die "Can't create 'inc_include_ref.txt' for writing. $!\n";
open my $non_include , '>', 'non_include.txt'
or die "Can't create 'non_include.txt' for writing. $!\n";
open my $overlap , '>', 'overlap.txt'
or die "Can't create 'overlap.txt' for writing. $!\n";
open my $to , '>', 'ref_sorted.txt' or die "Can't create for writing. $!\n";
while (my ($chr, $v) = each %ref) {
# sort region using start sites, if equal, using end sites
my @ref_data = sort {$a->[2] <=> $b->[2] or $a->[3] <=> $b->[3]} @$v;
# output sorted refGenes
print {$to} join("\t", @$_) . "\n" for @ref_data;
# if no such chromosome in inc, ignore and try next chr
next unless exists $data{$chr};
# sort region using start sites, if equal, using end sites
my @inc_data = sort {$a->[2] <=> $b->[2] or $a->[3] <=> $b->[3]} @{$data{$chr}};
my @store; # use to store some ref regions, since inc regions may overlapped
for (@inc_data) {
my @id = @$_;
if (@ref_data) {
my $result; # one of ref_include_inc, inc_include_ref, overlap or undef
while (@ref_data) {
my @rd = @{+shift @ref_data};
# igore genes on the left side of inc region
next if $rd[3] < $id[2];
# record each overlapped genes
push @store, [@rd];
# stop when this gene is on the right side of inc region
last if $rd[2] > $id[3];
# this gene must overlap with this inc region
if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
$result = 'ref_include_inc';
}
elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
$result = 'inc_include_ref';
}
else {
$result = 'overlap';
}
}
# output result
if ($result) {
if ($result eq 'ref_include_inc') {
print {$ref_include_inc} join("\t", @id), "\n";
}
elsif ($result eq 'inc_include_ref') {
print {$inc_include_ref} join("\t", @id), "\n";
}
else {
print {$overlap} join("\t", @id), "\n";
}
}
else { # if $result is undef, means no overlapped genes
print {$non_include} join("\t", @id), "\n";
}
# put this genes back into @ref_data, use to check next inc region
unshift @ref_data, @store;
@store = ();
}
else {
# no regions in ref left, print all regions left in inc to non_include
print {$non_include} join("\t", @id), "\n";
next;
}
}
delete $data{$chr};
}
while (my ($chr, $v) = each %data) {
my @inc_data = @$v;
for (@inc_data) {
print {$non_include} join("\t", @$_), "\n";
}
}
close $to;
close $ref_include_inc;
close $inc_include_ref;
close $non_include;
close $overlap;

复制代码

作者: yanglc2013 时间: 2013-03-28 20:43
回复 35# picbhan

这个程序思路确实很好，希望能加上注释。这样比较容易看。

作者: yanglc2013 时间: 2013-03-28 20:51
回复 35# picbhan

抱歉啊，刚复制到editplus才看到有注释，非常感谢啊。

作者: picbhan 时间: 2013-03-29 09:42
回复 36# yanglc2013

还可以更改一下，在判断是否属于包含还是overlap的时候加上last，这样可以在同一个inc的region与refgene存在多种情况时优先选择包含的情况，也可以减少运行时间。

if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
$result = 'ref_include_inc';
last;
}
elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
$result = 'inc_include_ref';
last;
}

复制代码

作者: yanglc2013 时间: 2013-03-29 14:49
回复 38# picbhan

大哥，输出结果中，有两种情况（non_include和overlap）少了数据，我没找到原因。请教你下。

test for help.rar (3.17 KB, 下载次数: 3)

作者: dahe_1984 时间: 2013-03-29 15:06

picbhan 发表于 2013-03-28 09:28
回复 8# dahe_1984

我自己写的是把ref中ch*作为Key值，然后foreach($ref[2]..$ref[3]),$hash{$[1]}++;

Inc中判断exists($hash{}), 然后判断in，out，还是overlap，生成的时候就跑死了，跟您的程度无关。

作者: picbhan 时间: 2013-03-29 18:28
本帖最后由 picbhan 于 2013-03-29 18:29 编辑

回复 40# dahe_1984

我只能说你太暴力，比inc中每一个region对ref中整个chr*的比较还要暴力。呵呵

作者: picbhan 时间: 2013-03-29 18:49
回复 39# yanglc2013

不是因为少了，而是inc中每个region我都确保只输出到一个文件中，没有输出到多个文件中，如果某个region与ref中的gene存在完全包含关系，就不会输出到overlap中。比如你说overlap中少了 lnc-TAF12-1:4 chr2 9 15 - ，那是因为ref中 NM_001097613 chr2 11 13 - 1697完全落在了inc的区间内，所以那个region就只输出到了inc_include_ref这个文件中。你如果要想一个region输出到多个文件中（对每种情况都输出）的话，可以稍微改一下代码就行，把标量$result改为%result,然后用不同情况下的字符串作为key，最近依据key值输出到不同文件中就行。但是这种情况下不能在找到包含关系之后就立即last（就是我最后帮你加的那两个last语句）。

另外你说non_include 少了 lnc-HES5-1:2 chr3 14 100 - ，这是因为你说不需要考虑正负链，那么这个region与ref中 NM_016505 chr3 50 128 + 827 是overlap的，所以输出在overlap中。
你再好好检查下吧。

作者: rubyish 时间: 2013-03-30 09:34
本帖最后由 rubyish 于 2013-03-30 05:35 编辑

ref

NM_000721 chr1 100 110 + 246
NM_020394 chr1 110 300 - 2470

lnc

lnc-FAM20C-5:2 chr1 100 150 +

输出文件1
OR
输出文件2

作者: yanglc2013 时间: 2013-03-30 09:54
回复 43# rubyish

lnc和ref比较，ref作为数据库。输出四个文件，输出的内容是lnc中的序列：

文件一：contain：就是ref图中的红色区间包含lnc图中红色区间。
文件二：overlap：就是lnc图中的红色区间和ref图中红色区间相交。
文件三：gapped：就是lnc图中的红色区间和ref图中红色区间包含，也不相交，是独立分开的。
文件四：reverse：就是lnc图中的红色区间包含ref图中的红色区间，和文件一刚好相反。

作者: yanglc2013 时间: 2013-03-30 10:07
回复 42# picbhan

大哥，其实是这样的，lnc中每个区间都要单独和ref中的区间比较，每一个染色体（比如lnc-HES5-1:2 chr3 14 100 -），它和NM_016505 chr3 50 128 + 827是overlap的，但是它和NR_026818 chr3 110 150 - 585却是non_include，所以输出的non_include文件中应该也包含lnc-HES5-1:2 chr3 14 100 - 啊。（就是lnc中同一个区间，他可能有很多种情况，都要分别输出。这也没有考虑正负链啊）

作者: yanglc2013 时间: 2013-03-30 10:22
回复 42# picbhan
我意思是对每种情况都输出啊。就是一个region，它可能是contain，可能是gapped，可能是overlap也可能lnc包含ref，如果这四种情况都存在，那么输出的四个文件中，都要包含这个region。大哥，真抱歉啊，我没表达清楚。你那个程序，我还没研究透，不是很会改你的程序。还是会出错。

作者: rubyish 时间: 2013-03-30 11:18
比如
lnc-HES5-1:2       chr3       14 100       -

NM_016505       chr3       50 128       +       827
NR_026818       chr3       110 150       -       585
NR_026835       chr3       12 150       -       586
NR_026822       chr3       16 50       -       534
NR_0268200    chr3       6 12       -       531

输出文件?

作者: yanglc2013 时间: 2013-03-30 14:10
本帖最后由 yanglc2013 于 2013-03-30 14:12 编辑

回复 47# rubyish
不是啊，输出的是lnc中的数据。（比如lnc-HES5-1:2 chr3 14 100 -）中的 [14,100]这个区间和ref中的（NM_016505 chr3 50 128 + 827
）中的[50,128]这个区间比较，它们有一个相交的区间[50,100]，所以在overlap文件中输出。如果和（NR_026818 chr3 110 150 - 585）比较，即区间[110,150]，没有公共部分，就输出到文件gapped。如果和（NR_026835 chr3 12 150 - 586）比较，lnc包含在ref区间中，输出到contain中，如果反过来，ref包含在lnc中，就输出到reserve文件中，输出的文件都是（lnc-HES5-1:2 chr3 14 100 -），如果四种情况都存在，那么，输出的四个文件中，都应该有（lnc-HES5-1:2 chr3 14 100 -）

作者: rubyish 时间: 2013-03-30 17:59
本帖最后由 rubyish 于 2013-03-31 12:53 编辑

{:3_188:}

Delete

有错误

复制代码

作者: rubyish 时间: 2013-03-31 07:41
本帖最后由 rubyish 于 2013-03-31 12:52 编辑

{:3_188:}

Delete

有错误

复制代码

作者: picbhan 时间: 2013-03-31 12:12
回复 45# yanglc2013

你觉得这样子输出有具体意义吗？inc中任何一个region都会被输出到non_include中去，只要你的refGene足够多。对于后面三种有交集的情况，如果要输出到多个文件中，你把下面的代码替换一下就行。

my %result; # one of ref_include_inc, inc_include_ref, overlap or undef
while (@ref_data) {
my @rd = @{+shift @ref_data};
# igore genes on the left side of inc region
next if $rd[3] < $id[2];
push @store, [@rd];
# stop when this gene is on the right side of inc region
last if $rd[2] > $id[3];
# this gene must overlap with this inc region
if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
$result{ref_include_inc} = 1;
}
elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
$result{inc_include_ref} = 1;
}
else {
$result{overlap} = 1;
}
}
# output result
if (%result) {
if ($result{ref_include_inc}) {
print {$ref_include_inc} join("\t", @id), "\n";
}
if ($result{inc_include_ref}) {
print {$inc_include_ref} join("\t", @id), "\n";
}
if ($result{overlap}) {
print {$overlap} join("\t", @id), "\n";
}
}
else { # if $result is undef, means no overlapped genes
print {$non_include} join("\t", @id), "\n";
}
# put this genes back into @ref_data, use to check next inc region
unshift @ref_data, @store;
@store = ();
}

复制代码

作者: yanglc2013 时间: 2013-04-01 18:59
回复 51# picbhan

不会出现这样情况的。因为ref中每个染色体，比如很多chr1都是来源于有一个DNA，这样发生的剪切，如果存在contain，那么另一种情况只能是gapped。

作者: dahe_1984 时间: 2013-04-02 15:57

rubyish 发表于 2013-03-31 07:41

你写的code怎么没有了？我还想学习下呢

作者: yanglc2013 时间: 2013-04-03 10:37
回复 51# picbhan

大哥，还有个问题请教你。
如果我想在输出的结果中的结果中，同时加上对应的ref中的数据，比如。
lnc-C1orf227-1:2 chr1 213028819 213031452 - NM_014053 chr1 213031596 213072705 + 2210

应该怎么修改print。

现在输出的结果形式如下：
lnc-C1orf227-1:2 chr1 213028819 213031452 -

作者: picbhan 时间: 2013-04-07 15:02
回复 54# yanglc2013

你在下面的代码中加入输出的信息（包括@rd和@id），并且把原来的三个赋值语句改为一个表示已经输出过的flag，这样在后面判断是否是non-overlap的时候就只需要看这个flag是不是真值了。

if ($id[2] >= $rd[2] && $id[3] <= $rd[3]) {
$result{ref_include_inc} = 1;
}
elsif ($rd[2] >= $id[2] && $rd[3] <= $id[3]) {
$result{inc_include_ref} = 1;
}
else {
$result{overlap} = 1;
}
if (%result) {
if ($result{ref_include_inc}) {
print {$ref_include_inc} join("\t", @id), "\n";
}
if ($result{inc_include_ref}) {
print {$inc_include_ref} join("\t", @id), "\n";
}
if ($result{overlap}) {
print {$overlap} join("\t", @id), "\n";
}
}
else { # if $result is undef, means no overlapped genes
print {$non_include} join("\t", @id), "\n";
}

复制代码

作者: yanglc2013 时间: 2013-04-08 19:17
回复 2# picbhan

谢谢啊，你的程序非常好用。还有个难题想请教你，我也写了几个，一塌糊涂啊。

难题.rar (1.89 MB, 下载次数: 9)

作者: picbhan 时间: 2013-04-09 13:46

yanglc2013 发表于 2013-04-08 19:17
回复 2# picbhan

谢谢啊，你的程序非常好用。还有个难题想请教你，我也写了几个，一塌糊涂啊。

首先，我没有看到你说的你也写了几个，我只看到你的问题。
其次，其实你的问题很简单，如果你真认真学了Perl并且自己动手写了，肯定是可以完成的，这根本不是什么难题，仅仅只是简单的格式转换而已。
最后，没有人可以一直帮你写代码，而且这种问题都等着别人给你写代码的话只会降低你自己的效率。好好学perl基础编程吧，足以在生物里用了。

#!perl
use List::MoreUtils 'pairwise';
$fin = 'refGene.txt';
open my $in, '<', $fin or die "Can't open $fin for reading. $!\n";
$fout = 'result.txt';
open my $to, '>', $fout or die "Can't create $fout for writing. $!\n";
while (<$in>) {
chomp;
my @t = split;
@start = split /,/, $t[1];
@end = split /,/, $t[2];
pairwise { print {$to} "$t[0]\t$a\t$b\n" } @start, @end;
}
close $in;
close $to;

复制代码

欢迎光临 Chinaunix (http://bbs.chinaunix.net/)