Project-UMIs/UMIsProject.R at master · mcmaniou/Project-UMIs · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
########## Loading libraries ##########
library(tidyverse)
library(data.table)
library(ShortRead)
library(stringr)
library(Biostrings)

source("casesWorkflows.R")
source("functions.R")

########## Inputs ##########

#type of data - paired or single
pairedData <- T

#UMI located in Read1 --> "R1"
#UMI located in Read1 and Read2 --> "R1 & R2"
UMIlocation <- "R1 & R2"

#length of the UMI
UMIlength <- 10

#length of th sequence
sequenceLength <- 251

#min read counts per UMI, for initial data cleaning
countsCutoff <- 0

#max UMI distance for UMI merging
UMIdistance <- 5

#max sequence distance for UMI correction
sequenceDistance <- 50

#inputs folder / working directory
inputsFolder <- "UMI in R1 and R2"

#outputs folder
outputsFolder <- "UMIc_output"

########## Run the appropriate scenario ##########


if (pairedData & UMIlocation == "R1"){   #case 1 -- paired data and UMI only in Read1

  inputFiles <- list.files(inputsFolder, pattern = "fastq")

  while(length(inputFiles)>0){

    file1 <- inputFiles[1]

    commonPart <- as.character(str_split(file1,"R1", simplify = T))
    commonPart <- commonPart[length(commonPart)]

    file2 <- str_replace(file1,paste0("R1",commonPart),paste0("R2",commonPart))

    filepath1 <- paste0(inputsFolder,"/",file1)
    filepath2 <- paste0(inputsFolder,"/",file2)

    pairedR1(filepath1, filepath2, outputsFolder,
             UMIlength, UMIdistance, sequenceLength, sequenceDistance, countsCutoff)


    inputFiles <- inputFiles[str_detect(inputFiles,paste0(file1,"|",file2), negate = T)]
  }


}else if (pairedData & UMIlocation == "R1 & R2"){   #case 2 -- paired data and UMI in Read1 and Read2

  inputFiles <- list.files(inputsFolder, pattern = "fastq")

  while(length(inputFiles)>0){

    file1 <- inputFiles[1]

    commonPart <- as.character(str_split(file1,"R1", simplify = T))
    commonPart <- commonPart[length(commonPart)]

    file2 <- str_replace(file1,paste0("R1",commonPart),paste0("R2",commonPart))

    filepath1 <- paste0(inputsFolder,"/",file1)
    filepath2 <- paste0(inputsFolder,"/",file2)

    pairedR1R2(filepath1, filepath2, outputsFolder,
               UMIlength, UMIdistance, sequenceLength, sequenceDistance, countsCutoff)


    inputFiles <- inputFiles[str_detect(inputFiles,paste0(file1,"|",file2), negate = T)]
  }


}else if (!pairedData){  #case 3 -- single data

  inputFiles <- list.files(inputsFolder, pattern = "fastq")

  while(length(inputFiles)>0){

    file1 <- inputFiles[1]

    filepath1 <- paste0(inputsFolder,"/",file1)

    single(filepath1, outputsFolder,
           UMIlength, UMIdistance, sequenceLength, sequenceDistance, countsCutoff)

    inputFiles <- inputFiles[str_detect(inputFiles,file1, negate = T)]
  }

}