Running mazter_mine with Single-end sequenced library

[ssun1116]

Dear Miguel,

Hello, Miguel. Recently I asked you a question about how to analyze a library sequenced by smRNA-Seq Kit.
When I asked you a question, at first I thought my library is suitable for paired-end analysis, but finally I realized that my library was sequenced in a way that is optimized only for single-end analysis.
My sequencing kit aims at analysis of diverse small RNA species, so the protocol strongly advices me that I should analyze only the read 1 in a single-end analysis.

You gave me advice that mapping only read1 will make problem of the end of fragment not mapped correctly. I fully agree with your advice and would like to follow your advice, but since my library has already been made in another lab, I have no choice but to use this data. In order to maintain the quality of the data as much as possible and to fully apply your analysis pipeline, I proceeded with the analysis by modifying your code little by little.

After many trials and errors, I was able to get a figure similar to the one in your paper, and I also got a cleavage efficiency table. But I'm a little worried whether my modification followed the analysis flow you created well.
Can I ask you to check if the modification I made was done correctly?

In the bam2ReadEnds.R code, I modified the samtools code :

## from
comm <- paste0("samtools view -bf 0x2 ", BAMfile," | bedtools bamtobed -i stdin -bedpe -mate1 | cut -f1-3,5-6,9 > ", BEDPfile)
## to
comm <- paste0("samtools view -bq 5 -F 4 ", BAMfile," | bedtools bamtobed -i stdin | cut -f1-3,6 > ", BEDPfile)

And, in the helperFunctions.R code, I added readBEDSE3 and singleEndReadsToGenes functions like this :

# Read BedSE file for getting read counts (strand specific)
readBEDSE3 <- function(fileName){
  tmp <- read.delim(fileName, header = F, stringsAsFactors = F)
  colnames(tmp) <- c("chr", "start", "end", "strand")
  tmp$start <- tmp$start + 1
  return(tmp)
}
# Function single end BEDfile to transcriptome
singleEndReadsToGenes <- function(parCluster, bedPEGzFile, gAnnot){
  inputFile <- readBEDSE3(gzfile(BEDPfile))
  commChrom <- intersect(unique(gAnnot$chr), unique(inputFile$chr))
  genPReads <- parLapply(parCluster, commChrom, function(iChr){
    tmpIF <- subset(inputFile, chr == iChr)
    tmpGA <- subset(gAnnot, chr == iChr)
    tmpO <- lapply(tmpGA$name, function(iGene){
      exBlock <- exonBlockGen(iGene, geneAnnot)
      if(geneAnnot[tmpGA$name, "strand"] == "+"){
          tmpDF <- subset(tmpIF, strand == "+" & start %in% exBlock & end %in% exBlock)
          tmpDF$start <- match(tmpDF$start, exBlock)
          tmpDF$end <- match(tmpDF$end, exBlock)
          tmpDF <- tmpDF[,c("chr", "start", "end")]
      }else if(geneAnnot[tmpGA$name, "strand"] == "-"){
          tmpDF <- subset(tmpIF, strand == "-" & start %in% exBlock & end %in% exBlock)
          tmpDF$start <- match(tmpDF$end, exBlock)
          tmpDF$end <- match(tmpDF$start, exBlock)
          tmpDF <- tmpDF[,c("chr", "start", "end")]
      }
      tmpDF <- ddply(tmpDF,.(chr, start, end), nrow)
      if(ncol(tmpDF) == 4){
        colnames(tmpDF)[4] <- "times"
        tmpDF$len <- tmpDF$end - tmpDF$start
      }
      rownames(tmpDF) <- NULL
      return(tmpDF)
    })
    names(tmpO) <- tmpGA$name
    return(tmpO)
  })
  # Flatten list
  tmp <- do.call(c, genPReads)
  genPReads <- tmp
  return(genPReads)
}

I think you'll be very busy reaching end of the PhD course as you mentioned before. I'm really sorry to post such a complicated and confusing question. I would be very grateful if you could take a look at this issue when you have some time.

Thank you for reading my post, and hope you have a good day.

Sincerly,
Seoyeon

[AngelCampos]

Hello Seoyeon, As you mention. I am currently in the last months of my PhD I find it difficult to be able to check this. Especially if I don't have access to the actual data. I recommend you corroborate your results by checking that in your data the known m6A sites are differentially methylated e.g. CleavageEfficiency of IME4-KO - ClvEff. of WT. Or at least that they have low cleavage efficiency compared to other similar sites. I will try to give it a check as soon as I have some time. Best

…

On Fri, Aug 6, 2021 at 11:50 AM Seoyeon Kim ***@***.***> wrote: Dear Miguel, Hello, Miguel. Recently I asked you a question about how to analyze a library sequenced by smRNA-Seq Kit. When I asked you a question, at first I thought my library is suitable for paired-end analysis, but finally I realized that my library was sequenced in a way that is optimized only for single-end analysis. My sequencing kit aims at analysis of diverse small RNA species, so the protocol strongly advices me that I should analyze only the read 1 in a single-end analysis. You gave me advice that mapping only read1 will make problem of the end of fragment not mapped correctly. I fully agree with your advice and would like to follow your advice, but since my library has already been made in another lab, I have no choice but to use this data. In order to maintain the quality of the data as much as possible and to fully apply your analysis pipeline, I proceeded with the analysis by modifying your code little by little. After many trials and errors, I was able to get a figure similar to the one in your paper, and I also got a cleavage efficiency table. But I'm a little worried whether my modification followed the analysis flow you created well. Can I ask you to check if the modification I made was done correctly? In the bam2ReadEnds.R code, I modified the samtools code : ## from comm <- paste0("samtools view -bf 0x2 ", BAMfile," | bedtools bamtobed -i stdin -bedpe -mate1 | cut -f1-3,5-6,9 > ", BEDPfile) ## to comm <- paste0("samtools view -bq 5 -F 4 ", BAMfile," | bedtools bamtobed -i stdin | cut -f1-3,6 > ", BEDPfile) And, in the helperFunctions.R code, I added readBEDSE3 and singleEndReadsToGenes functions like this : # Read BedSE file for getting read counts (strand specific) readBEDSE3 <- function(fileName){ tmp <- read.delim(fileName, header = F, stringsAsFactors = F) colnames(tmp) <- c("chr", "start", "end", "strand") tmp$start <- tmp$start + 1 return(tmp) } # Function paired end BEDfile to transcriptome singleEndReadsToGenes <- function(parCluster, bedPEGzFile, gAnnot){ inputFile <- readBEDSE3(gzfile(BEDPfile)) commChrom <- intersect(unique(gAnnot$chr), unique(inputFile$chr)) genPReads <- parLapply(parCluster, commChrom, function(iChr){ tmpIF <- subset(inputFile, chr == iChr) tmpGA <- subset(gAnnot, chr == iChr) tmpO <- lapply(tmpGA$name, function(iGene){ exBlock <- exonBlockGen(iGene, geneAnnot) if(geneAnnot[tmpGA$name, "strand"] == "+"){ tmpDF <- subset(tmpIF, strand == "+" & start %in% exBlock & end %in% exBlock) tmpDF$start <- match(tmpDF$start, exBlock) tmpDF$end <- match(tmpDF$end, exBlock) tmpDF <- tmpDF[,c("chr", "start", "end")] }else if(geneAnnot[tmpGA$name, "strand"] == "-"){ tmpDF <- subset(tmpIF, strand == "-" & start %in% exBlock & end %in% exBlock) tmpDF$start <- match(tmpDF$end, exBlock) tmpDF$end <- match(tmpDF$start, exBlock) tmpDF <- tmpDF[,c("chr", "start", "end")] } tmpDF <- ddply(tmpDF,.(chr, start, end), nrow) if(ncol(tmpDF) == 4){ colnames(tmpDF)[4] <- "times" tmpDF$len <- tmpDF$end - tmpDF$start } rownames(tmpDF) <- NULL return(tmpDF) }) names(tmpO) <- tmpGA$name return(tmpO) }) # Flatten list tmp <- do.call(c, genPReads) genPReads <- tmp return(genPReads) } I think you'll be very busy reaching end of the PhD course as you mentioned before. I'm really sorry to post such a complicated and confusing question. I would be very grateful if you could take a look at this issue when you have some time. Thank you for reading my post, and hope you have a good day. Sincerly, Seoyeon — You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub <#7>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACHMOKLYYOE6FW5VRP6UIUTT3OO6LANCNFSM5BVPLN4Q> .

-- *Miguel Angel García Campos* Website <https://angelcampos.github.io/> | Twitter <https://twitter.com/GarciaCamposM> | GScholar <https://scholar.google.com/citations?user=-J7nzzcAAAAJ&hl> | Instagram <https://www.instagram.com/garciacampos_m/> PhD candidate at the SchwartzLab <https://www.weizmann.ac.il/molgen/Schwartz/>, Molecular Genetics Department, Weizmann Institute of Science.