Как вычислить корреляции между подматрицами? - PullRequest
Новая
       10

Как вычислить корреляции между подматрицами?

0 голосов
/ 27 марта 2019

Мне нужно вычислить корреляции между столбцами моей матрицы, которая имеет более 800 000 строк в R. Я решил разбить эту матрицу на подматрицы (с 60 000 строк для каждой) и вычислить попарные корреляции между этими подматрицами.

Я использую SLURM.Я хочу распределить вычисление корреляции между двумя подматрицами на узле кластера, который я использую, для распараллеливания.

На данный момент я создал аргумент, который учитывает количество столбцов основной матрицы, которую я хочу вычислить.

Например, с помощью командной строки data = data [1: opt $ subset,] (в R) я могу вычислить корреляцию от моего 1-го столбца до 10.000-го: для этого я настроилмассивы в моем коде SLurm: subset = $ ((SLURM_ARRAY_TASK_ID * 10000)).Я определил 10 массивов, и поэтому первый вычислит вычисления от 1-го до 1 * 10000-го столбца, второй от 1-го до 2 * 10000 = 20 000-го столбца ....

С этим аргументом data = data [(as.numeric (opt $ subset) -4999): opt $ subset,], я могу вычислить корреляции в блок / подматрицу с определенным количеством столбцов.Например, если я хочу создать блоки из 5000 столбцов, я устанавливаю свой аргумент, как указано выше, и в SLURM с моими массивами: subset = $ ((SLURM_ARRAY_TASK_ID * 5000)).Итак, мой первый блок будет соответствовать от (1 * 5000) -4999 = 1-й столбец до 1 * 5000 = 5000-й столбец, второй блок будет соответствовать от (2 * 5000) -4999 = 5001-й столбец до 2 *5000 = 10.000-й столбец ..

Моя проблема здесь: корреляции вычисляются в эти блоки независимо.Что я хочу сделать, это вычислить корреляцию между всеми этими блоками следующим образом (= попарная корреляция между всеми блоками): 

        [,1] [,2]
  [1,]    1    1
  [2,]    1    2
  [3,]    1    3
  [4,]    1    4
  [5,]    1    5
  [6,]    1    6
 ...

до блока 6 между блоком 6.

Любой совет?

Приветствия

КОД R 

#load packages 
library(compositions)
library(parallel)
library(doParallel)
library(optparse)

args <- commandArgs(trailingOnly = F)

# get options 

option_list = list(
        make_option(c("-s", "--subset"), type="character", default=NULL, help="Input file matrix ")
);

opt_parser= OptionParser(usage = "Usage: %prog -f [FILE]",option_list=option_list, description= "Description:")

opt = parse_args(opt_parser)

#main code

print('Set Up Cores')

cores<-32
options('mc.cores'=cores)
registerDoParallel(cores)

print('Load matrice')

data<-read.table("/home/vipailler/PROJET_M2/raw/truelength2.prok2.uniref2.rares.tsv", sep="\t", h=T, row.names=1)+1

##THIS IS MY ARGUMENT###

#data=data[(as.numeric(opt$subset)-4999):opt$subset,]
data=data[1:opt$subset,]


res <- foreach(i = seq_len(ncol(data)),
 .combine = rbind,
 .multicombine = TRUE,
 .inorder = FALSE,
 .packages = c('data.table', 'doParallel')) %dopar% {
 if((i%%1000)==0){
 print(i)}
 apply(data, 2, function(x) 1 - ((var(data[,i] - x)) / (var(data[,i]) + var(x))))
}

КОД SLURM 

#!/bin/bash
#SBATCH --nodes=1
#SBATCH -o slurmjob-%A-%a.out
#SBATCH --job-name=rho_blocks_5k
#SBATCH --mail-user vincentpailler@hotmail.fr
#SBATCH --partition=normal
#SBATCH --time=1-00:00:00
#SBATCH --mem=250G
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --array=0-10

echo tableau de jobs numero $SLURM_ARRAY_JOB_ID, indices de $SLURM_ARRAY_TASK_MIN à $SLURM_ARRAY_TASK_MAX

echo $SLURM_ARRAY_TASK_ID

#Set up whatever package we need to run with

module load gcc/8.1.0 openblas/0.3.3 R

# SET UP DIRECTORIES

OUTPUT="$HOME"/PROJET_M2/bin/propr/$(date +"%Y%m%d")_parallel_blocks_32cpus_5000
mkdir -p "$OUTPUT"

export FILENAME=/home/vipailler/PROJET_M2/bin/coefficient_rho.R

subset=$((SLURM_ARRAY_TASK_ID*10000))

#Run the program

echo "Start job :"`date` >> "$OUTPUT"/temp_"$SLURM_ARRAY_TASK_ID".txt
echo "Start job :"`date`

Rscript $FILENAME --subset $subset  > "$OUTPUT"/"$SLURM_ARRAY_TASK_ID"

echo "Stop job : "`date` >> "$OUTPUT"/temp_"$SLURM_ARRAY_TASK_ID".txt
echo "Stop job : "`date`

Вывод, который я получаю:вот это: 

OTU0001     OTU0004    OTU0014    OTU0016    OTU0017      OTU0027
OTU0001  1.00000000  0.96688301 0.80621218 0.16754758 0.40818524  0.155976198
OTU0004  0.96688301  1.00000000 0.81330915 0.18928670 0.43247749  0.187540302
OTU0014  0.80621218  0.81330915 1.00000000 0.23753965 0.57237416  0.222890740
OTU0016  0.16754758  0.18928670 0.23753965 1.00000000 0.64007329  0.775772234
OTU0017  0.40818524  0.43247749 0.57237416 0.64007329 1.00000000  0.445145905
OTU0027  0.15597620  0.18754030 0.22289074 0.77577223 0.44514590  1.000000000
...

После этого я переупорядочив вывод: 

Df<-data.frame(var1=rownames(res)[row(res)[upper.tri(res)]],
        var2=colnames(res)[col(res)[upper.tri(res)]],
        corr=res[upper.tri(res)])

, чтобы получить: 

       var1    var2          corr
1   OTU0001 OTU0004  0.9668830120
2   OTU0001 OTU0014  0.8062121821
3   OTU0004 OTU0014  0.8133091522
4   OTU0001 OTU0016  0.1675475819
5   OTU0004 OTU0016  0.1892866996
6   OTU0014 OTU0016  0.2375396470
7   OTU0001 OTU0017  0.4081852433
8   OTU0004 OTU0017  0.4324774863
...
...