Count number of common columns between every row of a dataframe to create an all-vs-all matrix

First, let's keep only the PRESENT genes and build a frozenset of present genes per row:

s = (
 df.where(df.eq('PRESENT'))
   .stack()
   .reset_index(level=1)
   .groupby(level=0)['level_1'].agg(frozenset)
)

# 1    (gene3, gene1)
# 2    (gene2, gene1)
# 3    (gene3, gene2)
# Name: level_1, dtype: object

Now we can get the itertools.product (or, to be faster the combinations) to build the expected dataframe from the length of the intersection of genes set:

from itertools import product

out = (
 pd.Series({(i,j): len(a.intersection(b))
           for (i,a),(j,b) in  product(zip(s.index, s), repeat=2)})
   .unstack()
)

#    1  2  3
# 1  2  1  1
# 2  1  2  1
# 3  1  1  2

combinations are faster as only 1 of A/B or B/A is computed, and also not A/A:

from itertools import combinations

out = (
 pd.Series({(i,j): len(a.intersection(b))
           for (i,a),(j,b) in combinations(zip(s.index, s), r=2)})
   .unstack()
)

#      2    3
# 1  1.0  1.0
# 2  NaN  1.0

So I came up with another solution as well if anyone wants to use it;

Instead of using a dataframe with 'PRESENT' and 'LOST' I have used sed to transfer them into 1's and 0's.

sed -i '' 's/PRESENT/1/g' pav-matrix.binary.txt
sed -i '' 's/LOST/0/g' pav-matrix.binary.txt

Then used this binary matrix to compare every product of each row's and get the sum of the product array.

for (rowInQuestion in 1:nrow(pavMatrix_transposeBinary)){
  for (rowBeingComparedTo in 1:nrow(pavMatrix_transposeBinary)){
    pavSelfSimilarity[rowInQuestion,rowBeingComparedTo] <- sum(pavMatrix_transposeBinary[rowInQuestion,]*pavMatrix_transposeBinary[rowBeingComparedTo,])
  }
}

Vectorizing will give a speed boost for larger datasets:

test_PAV <- as.data.frame(matrix(sample(c("PRESENT", "LOST"), 3e4, replace = TRUE), ncol = 150))

fCommonVec <- function(test_PAV) {
  m <- as.matrix(test_PAV) == "PRESENT"
  n <- nrow(m) - 1L
  test_self <- matrix(nrow = nrow(test_PAV), ncol = nrow(test_PAV))
  test_self[lower.tri(test_self)] <- rowSums(m[rep.int(1:n, n:1),] & m[sequence(n:1, 2:nrow(m)),])
  diag(test_self) <- rowSums(m)
  return(test_self)
}

Compare with the loop version:

fCommonLoop <- function(test_PAV) {
  pavMatrix_transposeBinary <- as.matrix(test_PAV) == "PRESENT"
  pavSelfSimilarity <- matrix(nrow = nrow(pavMatrix_transposeBinary), ncol = nrow(pavMatrix_transposeBinary))
  
  for (rowInQuestion in 1:nrow(pavMatrix_transposeBinary)){
    for (rowBeingComparedTo in 1:nrow(pavMatrix_transposeBinary)){
      pavSelfSimilarity[rowInQuestion,rowBeingComparedTo] <- sum(pavMatrix_transposeBinary[rowInQuestion,]*pavMatrix_transposeBinary[rowBeingComparedTo,])
    }
  }
  
  return(pavSelfSimilarity)
}

test_self <- fCommonVec(test_PAV)
pavSelfSimilarity <- fCommonLoop(test_PAV)

all.equal(test_self[lower.tri(test_self, diag = TRUE)], pavSelfSimilarity[lower.tri(pavSelfSimilarity, diag = TRUE)])
#> [1] TRUE

microbenchmark::microbenchmark(fCommonVec = fCommonVec(test_PAV),
                               fCommonLoop = fCommonLoop(test_PAV))
#> Unit: milliseconds
#>         expr      min       lq      mean   median       uq      max neval
#>   fCommonVec  25.5148  30.6751  36.76842  32.6326  36.0579  76.0979   100
#>  fCommonLoop 215.4542 226.2288 235.84177 229.6783 235.9802 304.1668   100

UPDATE

For really large datasets, here is a function that iterates by column to cut down on memory usage:

test_PAV <- as.data.frame(matrix(sample(c("PRESENT", "LOST"), 267*113000, replace = TRUE), nrow = 267))

fCommonVec2 <- function(test_PAV) {
  m <- as.matrix(test_PAV) == "PRESENT"
  idx1 <- (nrow(m) - 1L):1
  idx2 <- rep.int(1:(nrow(m) - 1L), idx1)
  idx3 <- sequence(idx1, 2:nrow(m))
  v <- integer(length(idx2))
  
  for (i in 1:ncol(m)) {
    v <- v + (m[idx2, i] & m[idx3, i])
  }
  
  test_self <- matrix(nrow = nrow(test_PAV), ncol = nrow(test_PAV))
  test_self[lower.tri(test_self)] <- v
  diag(test_self) <- rowSums(m)
  return(test_self)
}

It took about 40 seconds on my machine and a small fraction of my 8 GB of RAM.

system.time(test_self2 <- fCommonVec2(test_PAV))
#>    user  system elapsed 
#>   36.52    0.24   36.76

This for loop with awk should work.

for i in {2..4}; do 
   for j in {2..4};do 
      tmp=$(awk '{print $'$i',$'$j'}' input |sed 1d|awk '$1==$2 {print}'|grep -c "PRESENT");echo -en "$tmp\t"
   done 
   echo -en "\n" 
done >> test_self

Output:

2       1       1
1       2       1
1       1       2

You might have to add the headers by yourself, the order should be the same.

Thanks!!