搜索
查看: 3791|回复: 3

多个差异分析结果直接量量取交集并集

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-4-3 14:38:52 | 显示全部楼层 |阅读模式
这个问题来自于楼上学生咨询我的Jaccard 指数的问题,大家请自行搜索,我这里不介绍背景了。
直接上题目:
文件如下,是我用perl单行命令模拟的数据:
[Perl] 纯文本查看 复制代码
perl -e 'BEGIN{use List::Util qw/shuffle/;@gene=A..Z}{foreach(1..10){@this_genes=@gene[(shuffle(0..25))[0..int(rand(20))+4]];print "comparison_$_ <-  ";print join(";",@this_genes);print "\n"}}'

解释起来也很简单,就是总共10次差异分析,每次差异分析得到的差异基因在同一行的后面用大小字母表示。

[AppleScript] 纯文本查看 复制代码
comparison_1 -> I;G;E;V;C;K;B;W
comparison_2 -> G;E;N;H;Y;M;L;S;K;A;J;O;D;P;R;U;Q;F;Z;C
comparison_3 -> Y;V;U;N;H;K;I;P;S;F;D;X;G;C;Z;J;Q;T;W;O;E;M
comparison_4 -> N;T;K;B;H;Z;W;C;Q;I;V;F;D;S;R;Y;J;X;P;O;G;L;A
comparison_5 -> G;J;A;H;W;T;Z;E;Y;S;R
comparison_6 -> Z;M;D;R;P;G;L;W;Y;U;X;E;A;S;T;I;H
comparison_7 -> H;Z;T;O;W;Q;M;X;J;N;U;K;F;P;I;C;S;Y;A;B
comparison_8 -> A;R;L;T;W;Q;S;F;P;X;E;V;Y;G;K;J;Z;C
comparison_9 -> J;X;K;D;O;H;L;F;C;P;R;N
comparison_10 -> G;S;K;H;C;O;W;F;Q;X


那么,作业就是求出每两个差异分析的结果的基因的交集的个数除以它们基因的并集的个数的比值。
得到一个比值矩阵。
我把我的R代码放在下面,希望你们用perl或者python来做出来。
[AppleScript] 纯文本查看 复制代码
a=readLines('test.txt')
n=unlist(lapply(a , function(x){
  trimws(strsplit(x,'<-')[[1]][1])
}))
v=lapply(a , function(x){
  trimws(strsplit(trimws(strsplit(x,'<-')[[1]][2]),';')[[1]])
})
names(v)=n
tmp=unlist(lapply(v, function(x){
  lapply(v, function(y){
    p1=length(intersect(x,y)) 
    p2=length(union(x,y)) 
    return(p1/p2)
  })
}))
out=matrix(tmp,nrow = length(n))
rownames(out)=n
colnames(out)=n
out[1:5,1:5]
结果的前五行如下:
[AppleScript] 纯文本查看 复制代码
             comparison_1 comparison_2 comparison_3 comparison_4 comparison_5
comparison_1    1.0000000    0.1666667    0.3043478    0.2916667    0.1875000
comparison_2    0.1666667    1.0000000    0.6800000    0.6538462    0.4090909
comparison_3    0.3043478    0.6800000    1.0000000    0.7307692    0.3750000
comparison_4    0.2916667    0.6538462    0.7307692    1.0000000    0.4166667
comparison_5    0.1875000    0.4090909    0.3750000    0.4166667    1.0000000





上一篇:用RSAT网页版工具在线寻找motif
下一篇:使用MEM数据库来做长链非编码RNA的靶基因预测
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-4-4 01:49:00 | 显示全部楼层
perl代码:

[Perl] 纯文本查看 复制代码
#!/usr/bin/perl -w
use strict;

my $in = "test.txt";
open my $fh0, $in or die "Cannot open $in";
while(<$fh0>){
        chomp;
        my @array = split /\s/,$_;
        my $header = $array[0];
        print "\t$header";
}
print "\n";
close $fh0;

open my $fh, $in or die "Cannot open $in";
while(<$fh>){
        chomp;
        my @array = split /\s/,$_;
        my $header = $array[0];
        print "$header\t";
        my @gene = split /;/, $array[2];
        open my $fh2, $in or die "Cannot open $in";
        while(<$fh2>){
                chomp;
                my %hash1=();
                my %hash2=();
                my @tmp1 = split /\s/, $_;
                my @tmp2 = split /;/, $tmp1[2];
                my @gene_all = (@tmp2,@gene);
                my @inter = grep {
                        ++$hash1{$_} > 1;
                }@gene_all;
                my @union = grep {
                        ++$hash2{$_} < 2;
                }@gene_all;
                my $num = sprintf "%.5f", (scalar @inter)/(scalar @union);
                print "$num\t";
        }
        print "\n";
        close $fh2;
}
close $fh;
回复 支持 反对

使用道具 举报

0

主题

2

帖子

123

积分

注册会员

Rank: 2

积分
123
发表于 2017-4-18 16:45:46 | 显示全部楼层
python 上三角形式
[Python] 纯文本查看 复制代码
#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""
test.py
============
Author:
    Chunyu

Description (zh-cn):
    Context
    
Reurirements:
    Python packages: 
"""
from __future__ import division
import pandas as pd
import numpy as np
def main():
    """ main """
    listval=[]
    listname=[]
    with open('diff_analysis','r') as input_file:
        for line in input_file.readlines():
            line=line.strip().split('\t')
            listname.append(line[0])
            listval.append(line[1].split(';'))
    arr=np.array(listval)
    rows=len(arr)
    tmp_array=np.zeros((rows,rows))
    for i in xrange(rows):
        for j in xrange(i,rows):
            tmp_array[i,j]=require_fun(arr[i],arr[j])
    df=pd.DataFrame(tmp_array,columns=listname,index=listname)
    print df
    df.to_csv('tmp2.csv')
def require_fun(i,j):
    union_set=list(set(i).union(set(j)))
    intersection_set=list(set(i).intersection(set(j)))
    print len(intersection_set)/len(union_set)*1.0
    return len(intersection_set)/len(union_set)
 
if __name__ == "__main__":
    main()

回复 支持 反对

使用道具 举报

1

主题

16

帖子

124

积分

注册会员

Rank: 2

积分
124
发表于 2017-9-18 11:29:33 | 显示全部楼层
[Perl] 纯文本查看 复制代码
use List::Compare;
open F,"ref.txt";
for (1..5){
    $tem=<F>;
    $tem=~/(\w;.*)/;
    @tem=(split ";",$1);
    $hash{$_}=[@tem];
}
close F;
for $first(1..5){
    @temp1=@{$hash{$first}};
    for (1..5){
        @temp2=values $hash{$_};
        $lc=List::Compare->new(\@temp1,\@temp2);
        $jj_num=$lc->get_intersection();
        $bj_num=$lc->get_union();
        $result=sprintf "%0.7f",$jj_num/$bj_num;
        push @{$_},$result;
    }
}
print "\tcomparison_1\tcomparison_2\tcomparison_3\tcomparison_4\tcomparison_5\n";
for (1..5){
    $a=join "\t",@{$_};
    print "comparison_$_\t".$a."\n";
}


用了List::Compare包取交集并集函数,毕竟现成的,不用白不用
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-10-19 01:27 , Processed in 0.031252 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.