atac-seq

# gzip fragment gz file
for i in `awk '{print $1}' df_meta_fragment.bed`
do
cd $i
gzip -d *.gz 
aws s3 cp ./ s3://scrna-seq-database/scATAC-seq/${i} --recursive
rm *
cd ../
done

## make mongodb

for i in `awk '{print $1}' df_meta_adult_fragment.bed`
do
cd ${i}
sudo service mongod restart && \
aws s3 cp s3://scrna-seq-database/scATAC-seq/${i}/ ./ --recursive && \
sort -k 1,1 -k2,2n -S 75% --parallel=16 `ls | grep "_fragment"` > sorted_fragments.bed && sed -i 1i'chr\tstart\tend\tbarcode\tcounts\tstrand' sorted_fragments.bed && \
sed -i '/^#/d' sorted_fragments.bed && \
mongoimport --db test --collection ${i} --type tsv --file sorted_fragments.bed --fields chr,start,end,barcode,counts && \
rm * && \
cd ../
done

#sort -k4,4 -S 75% --parallel=16 sorted_fragments.bed > sorted_bybarcodes_fragments.bed
#awk '{print $4}' sorted_bybarcodes_fragments.bed | uniq -c > cell_id_metadata.txt  

python3 ~/scATAC_atlas/ref_coverage_v2.py pbmc_10x


## make coverage
for i in `awk '{print $1}' df_meta_adult_fragment.bed`
do
cd ${i}
sudo service mongod restart && \
aws s3 cp s3://scrna-seq-database/scATAC-seq/${i}/cell_id_metadata.txt ./ && \
aws s3 cp s3://scrna-seq-database/scATAC-seq/${i}/metadata.txt ./ && \
python3 ~/scATAC_atlas/ref_coverage.py $i && \
aws s3 cp ./coverage.npz s3://scrna-seq-database/scATAC-seq/${i}/ && \
rm * && \
cd ../
done

for i in `awk '{print $1}' df_meta_adult_fragment.bed`; do cd ${i}; sudo service mongod restart && aws s3 cp s3://scrna-seq-database/scATAC-seq/${i}/ ./ --recursive && sort -k 1,1 -k2,2n `ls | grep "_fragment"` > sorted_fragments.bed && sed -i 1i'chr\tstart\tend\tbarcode\tcounts\tstrand' sorted_fragments.bed && mongoimport --db test --collection ${i} --type tsv --file sorted_fragments.bed --fields chr,start,end,barcode,counts && python3 /home/ubuntu/scATAC_atlas/ref_coverage.py $i && aws s3 cp ./coverage.npz s3://scrna-seq-database/scATAC-seq/${i}/ && rm * && cd ../; done

for i in `awk '{print $1}' df_meta_adult_fragment.bed`
do
cd ${i}
sudo service mongod restart && \
aws s3 cp s3://scrna-seq-database/scATAC-seq/${i}/ ./ --recursive && \
sort -k 1,1 -k2,2n `ls | grep "_fragment"` > sorted_fragments.bed && sed -i 1i'chr\tstart\tend\tbarcode\tcounts\tstrand' sorted_fragments.bed && \
mongoimport --db test --collection ${i} --type tsv --file sorted_fragments.bed --fields chr,start,end,barcode,counts && \
python3 ~/scATAC_seq/ref_coverage.py $i && \
aws s3 cp ./coverage.npz s3://scrna-seq-database/scATAC-seq/${i}/ && \
rm * && \
cd ../
done

# generate chr sorted fragments file
cd ${i}
bedtools sort -i *_fragments.bed > sorted_fragments.bed

# make mongodb database from fragment file
sed -i 1i'chr\tstart\tend\tbarcode\tcounts\tstrand' sorted_fragments.bed # add column names to tsv file

# build mongodb datasets
mongoimport --db test --collection adipose_omentum_SM-ADYHB_1 (${i}) --type tsv --file sorted_fragments.bed --fields chr,start,end,barcode,counts

snaptools snap-pre --input-file=sorted_bybarcodes_fragments.bed --output-snap=pbmc_10k.snap --genome-name=hg38 --genome-size=hg38.chrom.sizes --min-mapq=10 --min-flen=0 --max-flen=1000 --keep-chrm=TRUE --keep-single=FALSE -keep-secondary=False --overwrite=True --max-num=1000000 --min-cov=500 --verbose=True
snaptools snap-add-bmat --snap-file=pbmc_10k.snap --bin-size-list 5000 --verbose=True
snaptools snap-add-pmat --snap-file=pbmc_10k.snap --peak-file=pbmc_granulocyte_sorted_10k_atac_peaks.bed --verbose=True

library(SnapATAC)
library(GenomicRanges)

chrom <- paste0("chr", c(1:22, "X", "Y"))
pbmc_10k <- createSnap("/home/ubuntu/scATAC_atlas/10x_genomics_same_cells/pbmc/pbmc_10k.snap", 
                    sample = "lung1", num.cores = 26)
barcodes <- read.table("/home/ubuntu/scATAC_atlas/10x_genomics_same_cells/pbmc/cell_id_metadata.txt",
                           stringsAsFactors = F)$V1
pbmc_10k <- pbmc_10k[which(pbmc_10k@barcode %in% barcodes),]
pbmc_10k <- addBmatToSnap(pbmc_10k, bin.size=5000, num.cores=26)
pbmc_10k <- makeBinary(pbmc_10k, mat="bmat")

peaks <- read.table("~/scATAC_atlas/10x_genomics_same_cells/pbmc/pbmc_granulocyte_sorted_10k_atac_peaks.bed",
                    stringsAsFactors = F)
peaks <- peaks[peaks$V1 %in% chrom, ]
write.table(peaks, "~/scATAC_atlas/10x_genomics_same_cells/pbmc/peaks.bed", row.names = F, col.names = F,sep = "\t", quote = F)

idy = queryHits(findOverlaps(pbmc_10k@feature, black_list.gr))
if(length(idy) > 0){pbmc_10k = pbmc_10k[,-idy, mat="bmat"]}
chr.exclude = seqlevels(pbmc_10k@feature)[grep("random|chrM", seqlevels(pbmc_10k@feature))]
idy = grep(paste(chr.exclude, collapse="|"), pbmc_10k@feature)
if(length(idy) > 0){pbmc_10k = pbmc_10k[,-idy, mat="bmat"]}

bin.cov <- Matrix::colSums(pbmc_10k@bmat)
bin.cutoff <- 11764-118
bin.mincutoff <- 118
idy <- which(bin.cov <= bin.cutoff & bin.cov >= bin.mincutoff)
pbmc_10k <- pbmc_10k[, idy, mat="bmat"]

bins.df = as.data.frame(pbmc_10k@feature)[,1:3]
bins.df <- bins.df[bins.df$seqnames %in% chrom, ]
write.table(bins.df,file = "scATAC_atlas/10x_genomics_same_cells/pbmc/bins.bed",append=FALSE,
            quote= FALSE,sep="\t", eol = "\n", na = "NA", dec = ".", 
            row.names = FALSE, col.names = FALSE, qmethod = c("escape", "double"),
            fileEncoding = "")

Matrix::writeMM(obj = pbmc_10k@bmat, file = "/home/ubuntu/scATAC_atlas/10x_genomics_same_cells/pbmc/pbmc_10k_bmat.mtx")






# Run rscripts
Rscript ref_atac_coverage.R $i
# index one of the column
m <- mongo(collection = args[1], 
           url= "mongodb://127.0.0.1:27017/?compressors=disabled&gssapiServiceName=mongodb")
m$index(add = '{"barcode" : 1}')

# query rows by cell barcode
getCoverage <- function(barcode_i){
  cell_tmp <- m$find(sprintf('{"barcode":"%s"}', barcode_i))
  cell_tmp <- cell_tmp[cell_tmp$chr %in% chrom, ]
  coverage <- bedTools.coverage(bed1 = chr1_annotation, bed2=cell_tmp)$V7
  return(coverage)
}

n.cores <- parallel::detectCores() - 6
my.cluster <- parallel::makeCluster(
  n.cores, 
  type = "PSOCK"
)
doParallel::registerDoParallel(cl = my.cluster)
length(barcodes)
coverage_list <- foreach(i = 1:50,  .combine = 'append') %dopar% {
  list(getCoverage(barcodes[i]))
}

parallel::stopCluster(cl = my.cluster)

saveRDS(coverage_list, file = "coverage.rds")

# upload coverage_list to s3
aws s3 cp ./coverage.rds s3://scrna-seq-database/scATAC-seq/${i}/
rm * 
cd ../

# liftover from hg19 to hg38
tmp <- makeGRangesFromDataFrame(dat, T)
chainObject <- import.chain("~/hg19ToHg38.over.chain")
tmp <- liftOver(tmp, chain = chainObject)

bedTools.sort<-function(functionstring="~/bedtools2/bin/sortBed",bed1,opt.string="")
{
  #create temp files
  a.file=tempfile()
  out   =tempfile()
  options(scipen =99) # not to use scientific notation when writing out
  
  #write bed formatted dataframes to tempfile
  write.table(bed1,file=a.file,quote=F,sep="\t",col.names=F,row.names=F)
  
  # create the command string and call the command using system()
  command=paste(functionstring,"-i",a.file,opt.string,">",out,sep=" ")
  cat(command,"\n")
  try(system(command))
  
  res=read.table(out,header=F)
  unlink(a.file);unlink(out)
  return(res)
}


bedTools.merge<-function(functionstring="~/bedtools2/bin/mergeBed",bed1,opt.string="")
{
  #create temp files
  a.file=tempfile()
  out   =tempfile()
  options(scipen =99) # not to use scientific notation when writing out
  
  #write bed formatted dataframes to tempfile
  write.table(bed1,file=a.file,quote=F,sep="\t",col.names=F,row.names=F)
  
  # create the command string and call the command using system()
  command=paste(functionstring,"-i",a.file,opt.string,">",out,sep=" ")
  cat(command,"\n")
  try(system(command))
  
  res=read.table(out,header=F)
  unlink(a.file);unlink(out)
  return(res)
}

bedTools.coverage<-function(functionstring="~/bedtools2/bin/coverageBed",bed1,bed2,opt.string="")
{
  #create temp files
  a.file=tempfile()
  b.file=tempfile()
  out   =tempfile()
  options(scipen =99) # not to use scientific notation when writing out
  
  #write bed formatted dataframes to tempfile
  write.table(bed1,file=a.file,quote=F,sep="\t",col.names=F,row.names=F)
  write.table(bed2,file=b.file,quote=F,sep="\t",col.names=F,row.names=F)
  
  # create the command string and call the command using system()
  command=paste(functionstring,"-a",a.file,"-b",b.file,opt.string,">",out,sep=" ")
  cat(command,"\n")
  try(system(command))
  
  res=read.table(out,header=F)
  unlink(a.file);unlink(b.file);unlink(out)
  return(res)
}


args <- commandArgs(trailingOnly = TRUE)
chrname <- paste0("chr", args[1])
chr1_fragment <- read.table(paste0("~/scRNA_seq/",chrname,".bed"), stringsAsFactors = F)
chr1_annotation <- read.table(paste0("~/scRNA_seq/",chrname,"_annotation.bed"), stringsAsFactors = F)
barcodes <- read.table("~/scRNA_seq/all_barcodes.txt", stringsAsFactors = F)$V1
coverage_list <- list()
for(i in 1:length(barcodes)){
  cell_tmp <- chr1_fragment[chr1_fragment$V4 == barcodes[i],]
  coverage_list[[i]] <- bedTools.coverage(bed1 = chr1_annotation, bed2=cell_tmp)$V7
}
saveRDS(coverage_list, file = paste0("~/scRNA_seq/",chrname,"_coverage.rds"))

for(ch in c(3:22, "X", "Y")){
  chrname <- paste0("chr", ch)
  edbx <- filter(EnsDb.Hsapiens.v86, filter = ~ seq_name == ch)
  all_gene_exon <- exonsBy(edbx, by="gene", columns=listColumns(EnsDb.Hsapiens.v86, "exon"), use.names=T)
  chr1_exon_bed <- list()
  for(i in 1:length(all_gene_exon)){
    if(sum(table(strand(all_gene_exon[[i]])) == 0) != 2){
      next
    }else{
      gene_list <- list(data.frame(chr = paste0("chr", seqnames(all_gene_exon[[i]])),
                                   start = start(all_gene_exon[[i]]),
                                   end = end(all_gene_exon[[i]]), strand = strand(all_gene_exon[[i]]),
                                   gene_name = all_gene_exon[[i]]$gene_name))
      chr1_exon_bed <- append(chr1_exon_bed, gene_list)
    }
  }
  chr1_exon_bed_sorted <- list()
  chr1_exon_bed_merged <- list()
  for(i in 1:length(chr1_exon_bed)){
    chr1_exon_bed_sorted[[i]] <- bedTools.sort(bed1 = chr1_exon_bed[[i]])
    chr1_exon_bed_merged[[i]] <- bedTools.merge(bed1 = chr1_exon_bed_sorted[[i]])
  }
  for(i in 1:length(chr1_exon_bed)){
    chr1_exon_bed_merged[[i]]$V4 <- chr1_exon_bed_sorted[[i]]$V4[1]
    chr1_exon_bed_merged[[i]]$V5 <- chr1_exon_bed_sorted[[i]]$V5[1]
  }
  for(i in 1:length(chr1_exon_bed_merged)){
    chr1_exon_bed_merged[[i]]$V6 <- "exon"
  }
  for(i in 1:length(chr1_exon_bed_merged)){
    colnames(chr1_exon_bed_merged[[i]]) <- c("chr", "start", "end", "strand", "gene_name", "annotation")
  }
  
  distance <- vector()
  for(i in 1:length(chr1_exon_bed_merged)){
    distance[i] <- chr1_exon_bed_merged[[i]]$end[nrow(chr1_exon_bed_merged[[i]])] - chr1_exon_bed_merged[[i]]$start[1]
  }
  
  sum(distance>1000000)
  
  chr1_exon_bed_merged <- chr1_exon_bed_merged[-which(distance>1000000)]
  
  chr1_intron_bed <- list()
  for(i in 1:length(chr1_exon_bed_merged)){
    if(nrow(chr1_exon_bed_merged[[i]]) == 1){
      next
    }else{
      df <- data.frame(
        chr = chrname,
        start = chr1_exon_bed_merged[[i]]$end[1:(nrow(chr1_exon_bed_merged[[i]])-1)] + 1,
        end = chr1_exon_bed_merged[[i]]$start[2:nrow(chr1_exon_bed_merged[[i]])] - 1,
        strand = chr1_exon_bed_merged[[i]]$strand[1],
        gene_name = chr1_exon_bed_merged[[i]]$gene_name[1],
        annotation = "intron"
      )
      chr1_intron_bed <- append(chr1_intron_bed, list(df))
    }
  }
  
  chr1_gene_bed <- list()
  for(i in 1:length(chr1_exon_bed_merged)){
    chr1_gene_bed[[i]] <- data.frame(chr = chrname, start = chr1_exon_bed_merged[[i]]$start[1], 
                                     end = chr1_exon_bed_merged[[i]]$end[nrow(chr1_exon_bed_merged[[i]])],
                                     strand = chr1_exon_bed_merged[[i]]$strand[1],
                                     gene_name = chr1_exon_bed_merged[[i]]$gene_name[1],
                                     annotation = "gene_body")
  }
  
  chr1_promoter_bed <- list()
  for(i in 1:length(chr1_gene_bed)){
    if(chr1_gene_bed[[i]]$strand == "+"){
      chr1_promoter_bed[[i]] <- data.frame(chr = chrname, start = chr1_gene_bed[[i]]$start-5000, 
                                           end = chr1_gene_bed[[i]]$start-1,
                                           strand = chr1_gene_bed[[i]]$strand[1],
                                           gene_name = chr1_gene_bed[[i]]$gene_name[1],
                                           annotation = "promoter")
    }else{
      chr1_promoter_bed[[i]] <- data.frame(chr = chrname, start = chr1_gene_bed[[i]]$end+1, 
                                           end = chr1_gene_bed[[i]]$end+5000,
                                           strand = chr1_gene_bed[[i]]$strand[1],
                                           gene_name = chr1_gene_bed[[i]]$gene_name[1],
                                           annotation = "promoter")
    }
  }
  
  chr1_gene_and_promoter_bed <- list()
  for(i in 1:length(chr1_gene_bed)){
    tmp <- c(chr1_promoter_bed[[i]]$start, chr1_promoter_bed[[i]]$end, chr1_gene_bed[[i]]$start, chr1_gene_bed[[i]]$end)
    chr1_gene_and_promoter_bed[[i]] <- data.frame(chr = chrname, start = min(tmp), 
                                                  end = max(tmp),
                                                  strand = chr1_gene_bed[[i]]$strand[1],
                                                  gene_name = chr1_gene_bed[[i]]$gene_name[1],
                                                  annotation = "gene_and_promoter")
  }
  
  chr1_gene_and_promoter_bed_combine <- do.call(rbind, chr1_gene_and_promoter_bed)
  chr1_gene_and_promoter_bed_combine <- bedTools.sort(bed1=chr1_gene_and_promoter_bed_combine)
  
  chr1_intergenic_bed <- data.frame(chr = chrname, start = c(1,chr1_gene_and_promoter_bed_combine$V3+1), 
                                    end = c(chr1_gene_and_promoter_bed_combine$V2-1, 248956422), strand = ".",
                                    gene_name = ".", annotation = "intergenic")
  
  final_chr1_exon <- do.call(rbind, chr1_exon_bed_merged)
  final_chr1_intron <- do.call(rbind, chr1_intron_bed)
  final_chr1_promoter <- do.call(rbind, chr1_promoter_bed)
  final_chr1_intergenic <- chr1_intergenic_bed
  
  chr1_annotation <- rbind(final_chr1_exon, final_chr1_intron, final_chr1_promoter, final_chr1_intergenic)
  chr1_annotation <- chr1_annotation[chr1_annotation$end > chr1_annotation$start, ]
  chr1_annotation <- bedTools.sort(bed1=chr1_annotation)
  write.table(chr1_annotation, paste0("~/scRNA_seq/",chrname,"_annotation.bed"), quote = F, row.names = F, col.names = F, sep = "\t")
}