生信编程15.多个差异分析结果直接取交集和并集
有一些五六年前的学生们都成长为了各个生物信息学相关公司的小领导，而且他们都有了自己的公众号，知乎号，也算是一番人物。最近他们跟我反馈面试找不到或者说很难直接考核筛选到认真干活的生信工程师，挺有意思的。让我想起来了早在生信技能树论坛创立之初我为了引流，而规划的200个生信工程师面试题。值得继续分享：
该问题的背景为Jaccard问题，首先我们来了解下Jaccard指数。
Jaccard系数
Jaccard index，又称为Jaccard相似系数（Jaccard similarity coefficient）用于比较有限样本集之间的相似性和差异性， Jaccard系数值越大，样本相似度越高。
定义:给定两个集合A和B，Jaccard系数定义为A与B的交集的大小与A与B的并集大小的比值，定义如下：
简单了解之后，进入 正题。
总共进行十次差异分析，得到每次差异分析的结果，然后求出每两个差异分析的结果的基因的交集个数除以它们的基因的并集的格式的比值，得到一个比值矩阵。
R语言实现
第一步，创建一个模拟数据，假设已经进行了十次差异分析
my_list <- list()
nums <- c(10,8,7,11,6,9)
for (i in seq(1:10)){
my_list[[i]]<-sample(LETTERS[seq(1:26)],sample(nums,1))
}
names(my_list) <- paste('comparison',c(1:10),sep = '_')
结果如下：
$comparison_1
[1] "C" "H" "A" "D" "J" "X" "V" "F" "R" "S"
$comparison_2
[1] "F" "A" "E" "G" "K" "O"
$comparison_3
[1] "S" "G" "R" "P" "E" "X" "Z" "Y" "C" "T" "J"
$comparison_4
[1] "W" "M" "D" "X" "R" "S" "U"
$comparison_5
[1] "I" "R" "J" "U" "O" "G" "N"
$comparison_6
[1] "D" "Z" "R" "S" "H" "B" "T" "X" "P" "K"
第二步，计算Jaccard系数
comp<-unlist(lapply(my_list, function(x){
lapply(my_list, function(y){
p1 = length(intersect(x,y))
p2 = length(union(x,y))
return(p1/p2)
})
}))
#将最终的结果转换为矩阵
comp_mat = matrix(comp, nrow = 10)
rownames(comp_mat) <- names(my_list)
colnames(comp_mat) <- names(my_list)
comp_mat[1:4,1:4]
结果文件如下：
comparison_1 comparison_2 comparison_3 comparison_4
comparison_1 1.0000000 0.1428571 0.3125000 0.3076923
comparison_2 0.1428571 1.0000000 0.1333333 0.0000000
comparison_3 0.3125000 0.1333333 1.0000000 0.2000000
comparison_4 0.3076923 0.0000000 0.2000000 1.00000000
R语言实现出来还蛮简单的，那我们再来试试python吧
Python实现
同样也分为两部分，第一部分创建模拟数据
#!/usr/bin/python
import pandas as pd
import numpy as np
import random
import string
#create data
my_list = {}
nums = list(range(6,11))
for i in list(range(1,11)):
letter = string.ascii_uppercase
#print(letter)
num = random.sample(nums,1)[0]
#print(num)
my_list['comparison_'+str(i)]=random.sample(letter,num)
结果文件如下：
comparison_1 ['C', 'P', 'U', 'X', 'O', 'E', 'L', 'Q', 'K']
comparison_2 ['E', 'V', 'G', 'N', 'S', 'Y', 'W', 'B', 'U']
comparison_3 ['A', 'X', 'Y', 'M', 'R', 'F', 'O', 'S', 'D']
comparison_4 ['R', 'L', 'K', 'N', 'Q', 'B', 'H']
comparison_5 ['V', 'O', 'U', 'F', 'P', 'M', 'D', 'Z', 'K']
comparison_6 ['T', 'Y', 'K', 'E', 'C', 'D', 'X', 'I', 'M']
comparison_7 ['M', 'E', 'U', 'W', 'O', 'R']
comparison_8 ['I', 'W', 'N', 'F', 'R', 'Y', 'Q', 'T']
comparison_9 ['I', 'Q', 'U', 'T', 'W', 'P']
comparison_10 ['D', 'C', 'Q', 'I', 'R', 'Y']
第二部分计算Jaccard系数
comp_list = []
for k1,v1 in my_list.items():
print(k1,v1)
for k2,v2 in my_list.items():
comp_name = k1+'vs'+k2
p1 = len(list(set(v1) & set(v2)))
p2 = len(list(set(v1) | set(v2)))
comp = round(p1/p2,3)
comp_list.append(comp)
mat= np.matrix(comp_list).reshape(10,10)
df = pd.DataFrame(mat)
df.columns=my_list.keys()
df.index = my_list.keys()
print(df)
结果文件如下：
comparison_1 comparison_2 ... comparison_9 comparison_10
comparison_1 1.000 0.125 ... 0.250 0.154
comparison_2 0.125 1.000 ... 0.154 0.071
comparison_3 0.125 0.125 ... 0.000 0.250
comparison_4 0.231 0.143 ... 0.083 0.182
comparison_5 0.286 0.125 ... 0.154 0.071
comparison_6 0.286 0.125 ... 0.154 0.364
comparison_7 0.250 0.250 ... 0.200 0.091
comparison_8 0.062 0.214 ... 0.400 0.400
comparison_9 0.250 0.154 ... 1.000 0.200
comparison_10 0.154 0.071 ... 0.200 1.000