Вычисление перекрытия генов между различными списками генов как% - PullRequest
0 голосов
/ 04 декабря 2018

Я сгенерировал таблицу, чтобы показать совпадение между различными списками генов.Так как у меня есть восемь различных списков генов, у меня есть 64 результата.Код, который у меня сейчас есть, выглядит следующим образом:

#-------------------------------------------------------------------------------
# Set the working directory and load the data files
#-------------------------------------------------------------------------------
setwd("~/Desktop/R_Project/Gene_overlap")
getwd()
files <- list.files(pattern="*.txt", full.names = TRUE)
files

data.list <- lapply(files, function(fil) {
  scan(file=fil, what=character())
})

names(data.list) <- basename(files) %>% stringr::str_remove("\\.txt$")

str(data.list)
# List of 8
# $ GSE108363_BCGdown_D: chr [1:350] "IL1B" "IL6" "IL1A" "CCL20" ...
# $ GSE108363_BCGdown_V: chr [1:267] "IL6" "CCL20" "IL1A" "CXCL5" ...
# $ GSE108363_BCGup_D  : chr [1:250] "FABP4" "CMTM2" "FUCA1" "CD36" ...
# $ GSE108363_BCGup_V  : chr [1:429] "FCN1" "FCGR3B" "MNDA" "CPVL" ...
# $ GSE108363_MTBdown_D: chr [1:86] "CCL20" "IL1B" "IL1A" "IL6" ...
# $ GSE108363_MTBdown_V: chr [1:244] "IL1B" "IL1A" "CCL20" "IL6" ...
# $ GSE108363_MTBup_D  : chr [1:128] "FUCA1" "FGL2" "TGFBI" "CPVL" ...
# $ GSE108363_MTBup_V  : chr [1:286] "FABP4" "RNASE1" "MNDA" "CPVL" ...

intersect(data.list$GSE108363_BCGdown_D, data.list$GSE108363_BCGdown_V) %>% length

sapply(data.list, length)



#-------------------------------------------------------------------------------
# Using the intersect function to see the overlaps 
#-------------------------------------------------------------------------------

data.file1 <- "GSE108363_BCGdown_V.txt"
data.file2 <- "GSE108363_BCGdown_D.txt"
data.file3 <- "GSE108363_BCGup_V.txt"
data.file4 <- "GSE108363_BCGup_D.txt"
data.file5 <- "GSE108363_MTBdown_V.txt"
data.file6 <- "GSE108363_MTBdown_D.txt"
data.file7 <- "GSE108363_MTBup_V.txt"
data.file8 <- "GSE108363_MTBup_D.txt"

genevect1 <- scan(data.file1, what=character(), sep="\n")
genevect2 <- scan(data.file2, what=character(), sep="\n")
genevect3 <- scan(data.file3, what=character(), sep="\n")
genevect4 <- scan(data.file4, what=character(), sep="\n")
genevect5 <- scan(data.file5, what=character(), sep="\n")
genevect6 <- scan(data.file6, what=character(), sep="\n")
genevect7 <- scan(data.file7, what=character(), sep="\n")
genevect8 <- scan(data.file8, what=character(), sep="\n")


filelist <- list(data.file1, data.file2, data.file3, data.file4, data.file5, data.file6, data.file7, data.file8)
all(sapply(filelist, file.exists))

# read files:
gene.lists <- lapply(filelist, function(f) {
  scan(file=f, what=character())
})


# set up empty matrix
x <- (length(gene.lists))^2
x
y <- rep(NA, x)
mx <- matrix(y, ncol=length(gene.lists))
mx
row.names(mx) <- sapply(filelist, basename) %>% stringr::str_remove('.txt$')
colnames(mx) <- sapply(filelist, basename) %>% stringr::str_remove('.txt$')
mx

mx.overlap.count <- mx

# seq_along(gene.lists) # 1 2 3 4 5 6 7 8
for (i in seq_along(gene.lists)) {
  g1 <- gene.lists[[i]]
  for (j in seq_along(gene.lists)) {
    g2 <- gene.lists[[j]]
    a <- intersect(g1, g2)
    b <- length(a)
    mx.overlap.count[j,i] <- b
  }
}

mx.overlap.count
View(mx.overlap.count)

Теперь я хотел бы сделать что-то похожее, но вместо просмотра перекрытия в виде чисел, я хотел бы увидеть степень существования перекрытий между различными списками генов. в процентах .Каким-то образом мне нужно будет узнать, больше ли g1 или g2, чтобы разделить b на большее, а затем умножить на 100. Любые предложения будут с благодарностью приняты.

1 Ответ

0 голосов
/ 04 декабря 2018

с использованием образца письма, поскольку вы не предоставили список генов:

set.seed(1)
data.list <- lapply(sample(10:20), function(n)LETTERS[sample(1:26, n)])

overlaps <- sapply(data.list, function(g1) 
  sapply(data.list, function(g2)
  {round(length(intersect(g1, g2)) / length(g2) * 100)}))

overlaps
      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
 [1,]  100   50   67   75   42   92   58   92   67    33    92
 [2,]   46  100   62   69   54   77   62   69   69    54    62
 [3,]   53   53  100   60   33   73   60   73   80    33    60
 [4,]   53   53   53  100   47   71   53   76   53    29    82
 [5,]   45   64   45   73  100   91   64   82   36    45    73
 [6,]   61   56   61   67   56  100   56   89   56    33    72
 [7,]   50   57   64   64   50   71  100   86   71    50    64
 [8,]   55   45   55   65   45   80   60  100   60    40    80
 [9,]   50   56   75   56   25   62   62   75  100    38    69
[10,]   40   70   50   50   50   60   70   80   60   100    70
[11,]   58   42   47   74   42   68   47   84   58    37   100

(я использовал set.seed, чтобы вы могли воспроизвести пример).Он использует вложенный объект sapply для итерации по обоим спискам генов по отдельности, а затем вычисляет процент для каждой комбинации векторов гена путем деления длины пересечения на общую длину второго вектора гена.Если вы хотите разделить на длину самого длинного из 2 генных векторов, замените length(g2) на max(length(g1), length(g2))

Добро пожаловать на сайт PullRequest, где вы можете задавать вопросы и получать ответы от других членов сообщества.
...