Skip to content

Commit

Permalink
a new example in folder T8 is added.
Browse files Browse the repository at this point in the history
  • Loading branch information
goknurginer committed May 17, 2024
1 parent 187e0e7 commit e99ad06
Show file tree
Hide file tree
Showing 25 changed files with 620,902 additions and 25,606 deletions.
1 change: 1 addition & 0 deletions .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -7,3 +7,4 @@
fastqs/
bams/
fastq_files/
index/
77 changes: 77 additions & 0 deletions R/counting_script_T8.R
Original file line number Diff line number Diff line change
@@ -0,0 +1,77 @@
# counting T8 example ----
# dir.create("T8")
path = "/Users/giner.g/Documents/Github/CuReSPR/T8"
## mageck ----
# dir.create(paste0(path,"/mageck"))
setwd(paste0(path,"/mageck"))
library(reticulate) # to run python code in r
venv_path <- "/Users/giner.g/miniconda3/envs/mageckenv/bin/python"
use_python(venv_path, required = TRUE)
system("/Users/giner.g/miniconda3/envs/mageckenv/bin/mageck count -l ../guidelibrary/brunello.tsv -n ./T8_results --sample-label Control,Vehicle,T8 --fastq ../fastqs/T0-Control.fastq.gz ../fastqs/T8-Vehicle.fastq.gz ../fastqs/T8-APR-246.fastq.gz", intern = TRUE)

## RSubread ----
library(Biostrings)
library(Rsubread)
library(GenomicAlignments)
# dir.create(paste0(path,"/rsubread"))
setwd(paste0(path,"/rsubread"))
guides <- read.csv("../guidelibrary/brunello.csv",header=FALSE)
head(guides)
sgRNAs <- DNAStringSet(guides$V2)
names(sgRNAs) <- guides$V1
sgRNAs
# dir.create("./index")
writeXStringSet(sgRNAs, file = "./index/guides.fa")
buildindex("./index/guides", "./index/guides.fa", indexSplit = FALSE)
setwd(paste0(path,"/fastqs"))
# single-end
reads <- dir(full.names = TRUE, pattern = "*.fastq*")
counts <- list()
mapping_results <- list()
for (i in 1:length(reads)) {
mapping_results[[i]] <- align("../rsubread/index/guides",
readfile1 = reads[i],
output_file = gsub(".fastq", ".bam",
tools::file_path_sans_ext(reads[i])),
nthreads = 4,
unique = TRUE,
nBestLocations = 1,
type = "DNA",
TH1 = 1,
maxMismatches = 0,
indels = 0)
temp <- readGAlignments(gsub(".fastq", ".bam",
tools::file_path_sans_ext(reads[i])))
counts[[i]] <- data.frame(table(seqnames(temp[width(temp) == "20"])), row.names = "Var1")
}
my_counts <- do.call(cbind, counts)
colnames(my_counts) <- reads
write.table(my_counts, paste0(path,"/rsubread/counts.txt"))

# paired-end
reads1 <- dir(full.names = TRUE, pattern = "*1.fastq*")
reads2 <- dir(full.names = TRUE, pattern = "*2.fastq*")
counts <- list()
mapping_results <- list()
for (i in 1:length(reads)) {
mapping_results[[i]] <- align("../rsubread/index/guides",
readfile1=reads1[i],
readfile2=reads2[i],
output_file = gsub(".fastq.gz", ".bam",
tools::file_path_sans_ext(reads[i])),
nthreads = 4,
unique = TRUE,
nBestLocations = 1,
type = "DNA",
TH1 = 1,
maxMismatches = 0,
indels = 0)
temp <- readGAlignments(gsub(".fastq.gz", ".bam",
tools::file_path_sans_ext(reads[i])))
counts[[i]] <- data.frame(table(seqnames(temp[width(temp) == "20"])), row.names = "Var1")
}
my_counts <- do.call(cbind, counts)
colnames(my_counts) <- reads
write.table(my_counts, paste0(path,"/rsubread/counts.txt"))

# wehi-marek
Loading

0 comments on commit e99ad06

Please sign in to comment.