R/conshaplotypes.R
In aliafdz/QApckg: Quality assessment for Miseq data derived from viral sequencing
Defines functions
ConsHaplotypes
Documented in
ConsHaplotypes
#' @title Generate consensus haplotypes
#'
#' @description Computes the intersection of forward and reverse strand haplotypes and generates some report files.
#'
#' @details This function is designed to be used after the execution of \code{\link{demultiplexPrimer}} function
#' from the same package. After the generation of FASTA files containing forward and reverse strand reads
#' for the evaluated samples, \code{\link{ConsHaplotypes}} executes multiple alignment with \code{\link{muscle}} and
#' returns the consensus haplotypes using \code{\link{IntersectStrandHpls}}, that will be saved using
#' the helper function \code{\link{SaveHaplotypes}}.
#'
#'
#' @param trimfiles Vector including the paths of demultiplexed files by specific primer, with fna extension.
#' @param pm.res The list returned by \code{\link{demultiplexPrimer}}, including \code{fileTable}
#' and \code{poolTable} data frames.
#' @param thr Threshold to filter haplotypes at minimum abundance before multiple alignment.
#' @param min.seq.len Threshold to filter haplotypes at minimum length before intersection.
#' @param max.difs Maximum number of mismatches allowed in resulting consensus haplotypes with
#' respect to the dominant one.
#'
#' @return The function returns a \code{\link{data.frame}} object containing the intersection results
#' for each combination of patient and amplicon region, including the initial number of reads, filtered out reads
#' (for being below a given frequency threshold or unique to a single strand), overlapping frequency
#' between both strands and the common reads (in percentage and nº of reads).
#'
#' @return After execution, two FASTA files for each combination of sample and pool will be saved in a newly
#' generated MACH folder; the first includes multiple alignment between forward and reverse strand haplotypes,
#' and the second includes the forward and reverse strands intersected.
#' Additionaly, some report files will be generated in the reports folder:
#' \enumerate{
#' \item{\code{MA.Intersects-SummRprt.txt}: Includes the sumary results by reads number after abundance filter and
#' strand intersection.}
#' \item{\code{MA.Intersects.plots.pdf}: Includes different barplots for each sample representing the frequency of
#' forward, reverse and intersected strand haplotypes.}
#' \item{\code{IntersectBarplots.pdf}: Includes different barplots for all combinations of patient and pool,
#' representing the number of intersected and filtered out reads, the intersection yield and global yield.}}
#'
#' @note A new file named \code{muscle.log} containing \code{\link{muscle}} options will be generated and saved in a folder named "tmp".
#'
#' @import Biostrings
#' @import ape
#' @import QSutils
#' @import parallel
#' @importFrom muscle muscle
#' @importFrom RColorBrewer brewer.pal
#' @importFrom methods as
#' @importFrom methods is
#' @seealso \code{\link{muscle}}, \code{\link{IntersectStrandHpls}}, \code{\link{demultiplexPrimer}}, \code{\link{SaveHaplotypes}}
#' @examples
#' splitDir <- "./splits"
#' # Save the file names with complete path
#' splitfiles <- list.files(splitDir,recursive=TRUE,full.names=TRUE,include.dirs=TRUE)
#' # Get data
#' samples <- read.table("./data/samples.csv", sep="\t", header=T,
#' colClasses="character",stringsAsFactors=F)
#' mids <- read.table("./data/mids.csv", sep="\t", header=T,
#' stringsAsFactors=F)
#' # Apply previous function from QA analysis
#' pm.res <- demultiplexPrimer(splitfiles,samples,primers)
#' # Save the files generated by previous function
#' trimDir <- "./trim"
#' trimfiles <- list.files(trimDir,recursive=TRUE,full.names=TRUE,include.dirs=TRUE)
#' # Define necessary parameters
#' min.seq.len <- 150
#' thr <- 0.2
#' int.res <- ConsHaplotypes(trimfiles, pm.res, thr, min.seq.len)
#'
#' @export
#' @author Alicia Aranda
ConsHaplotypes <- function(trimfiles,pm.res,thr=0.2, min.seq.len=150,max.difs=250){
# Si la ruta on es troben els fitxers de la carpeta trim no està ben especificada, intenta buscar la
# ruta correcta a partir del directori de treball
# Si tot i així no troba fitxers, atura l'execució i mostra un missatge d'error
if(length(trimfiles)==0) {
trimfiles <- list.files(paste(getwd(),"/trim",sep=''))
if(length(trimfiles)==0) {
stop("Please indicate correct path for demultiplexed reads by primer.\n")
}
}
# Si la taula pm.res no s'ha inclòs o no existeix, retorna un missatge d'error
if(missing(pm.res)|!exists('pm.res')|length(pm.res)==0) {
stop("The list obtained by demultiplexPrimer function is needed.\n")
}
if(!is.list(pm.res)){
stop("pm.res argument must be a list.\n")
}
# Si cap dels fitxers indicats a la carpeta trim no existeix, atura l'execució i
# mostra un missatge d'error
if(any(!file.exists(trimfiles))){
stop(paste(trimfiles[!file.exists(trimfiles)],"does not exist.\n"))
}
# Si no existeixen les carpetes on es guarden els fitxers resultants de la funció,
# es generen automàticament a la carpeta de treball
if(!dir.exists("./MACH")) {
dir.create("./MACH")}
mach.Dir <- "./MACH"
if(!dir.exists("./reports")) {
dir.create("./reports")}
repDir <- "./reports"
if(!dir.exists("./tmp")) {
dir.create("./tmp")}
# En lloc de carregar el RData es carreguen les taules fileTable i poolTable de la variable de la funció anterior.
FlTbl<- pm.res$fileTable
PoolTbl<- pm.res$poolTable
# Retorna els indexs de la taula FlTbl derivats d'ordenar les mostres en funció de l'ID del
# pacient, de l'amplicó (regió) avaluat i la cadena forward o reverse
o <- order(FlTbl$Pat.ID,FlTbl$Ampl.Nm,FlTbl$Str)
# Reordena les entrades de la taula en funció dels indexs anteriors, de manera que s'agrupen
# les entrades corresponents al mateix pacient i en segon ordre a la mateixa regió avaluada
FlTbl <- FlTbl[o,]
### Indica quins dels fitxers FASTA inclosos al directori trim corresponen a
# cadascuna de les cadenes
idx.fw<-grep("PrFW",trimfiles)
idx.rv<-grep("PrRV",trimfiles)
# Guarda els noms dels fitxers resultants d'aquest script, que correspondran a la concatenació
# del terme MACHpl02, l'ID del pacient i el nom de la regió avaluada, indicats de manera que només
# s'obtingui un fitxer .fna per a cada regió del mateix pacient (per això el condicional fw)
out.flnms <- paste("MACHpl02",FlTbl$Pat.ID[idx.fw],
FlTbl$Ampl.Nm[idx.fw],"fna",sep=".")
# Guarda la ruta on es desaran els fitxers .fna, a la carpeta MACH
out.flnms <- file.path(mach.Dir,out.flnms)
# Guarda els noms d'uns altres fitxers resultants, que correspondran a la concatenació
# del terme MAfwrv, l'ID del pacient i el nom de la regió avaluada, indicats de manera que només
# s'obtingui un fitxer .fna per a cada regió del mateix pacient (per això el condicional fw)
ma.flnms <- paste("MAfwrv",FlTbl$Pat.ID[idx.fw],
FlTbl$Ampl.Nm[idx.fw],"fna",sep=".")
# Guarda la ruta on es desaran els fitxers .fna, a la carpeta MACH
ma.flnms <- file.path(mach.Dir,ma.flnms)
# Guarda la meitat de la longitud de fitxers presents a la carpeta trim del fitxer previ
n <- length(trimfiles)/2
# Genera un data frame inicialment buit, amb n files (2 entrades per pacient) i 5 columnes
# que inclouran els resultats obtinguts en pasos posteriors per a les cadenes fw
rdf.fw <- data.frame(fw.all=integer(n),fw.lowf=integer(n),
fw.in=integer(n),fw.unq=integer(n),fw.com=integer(n))
# Genera un altre data frame també buit, amb n files i 5 columnes per incloure resultats
# de les cadenes rv
rdf.rv <- data.frame(rv.all=integer(n),rv.lowf=integer(n),
rv.in=integer(n),rv.unq=integer(n),rv.com=integer(n))
# Genera un altre data frame buit, amb n files i 6 columnes, per incloure els resultats globals
rdf.gbl <- data.frame(all=integer(n),lowf=integer(n),unq=integer(n),
ovrlp=numeric(n),common=numeric(n),Fn.rd=integer(n))
# Genera el primer fitxer pdf on s'inclouran diverses representacions de resultats generats en el bucle for
# Per cada regió amplificada de cada pacient, es representen els haplotips de cada cadena i els que han
# coincidit entre les dues, amb les corresponents freqüències
pdf(file.path(repDir,"MA.Intersects.plots.pdf"),paper="a4",
width=7,height=11)
par(mfrow=c(3,1))
# Bucle per iterar sobre totes les mostres (2 per pacient)
foreach(i=c(1:n)) %do% {
# Aplica la funció 'ReadAmplSeqs()' del paquet QSutils que permet obtenir els haplotips i les
# seves freqüències a partir d'un fitxer FASTA indicat.
lstf <- ReadAmplSeqs(trimfiles[idx.fw[i]])
# Guarda les seqüències de longitud major al mínim permès i filtra el nº de reads
# eliminant els de les seqüències que presenten longitud menor al mínim permès
lstf$nr <- lstf$nr[nchar(lstf$hseqs) >= min.seq.len]
# Aplica la funció 'GetQSData()' del paquet QSutils per llegir els amplicons del
# fitxer amb els valores d'abundància, filtrar per mínima abundància i ordenar
# segons les mutacions.
s1 <- GetQSData(trimfiles[idx.fw[i]],thr)
# Guarda les seqüències de longitud major al mínim permès
ffilt <- nchar(s1$seqs) >= min.seq.len
# Filtra les seqüències per eliminar les que presenten longitud menor al mínim permès
seqsf <- s1$seqs[ffilt]
nrf <- s1$nr[ffilt]
# Calcula la freqüència relativa de cada amplicó filtrat per abundància respecte el total
# de reads del fitxer inicial.
frq_lst <- round(nrf/sum(lstf$nr)*100,2)
# Donat que la funció del paquet QSutils retorna els noms dels haplotips diferent al desitjat, es canvien els
# noms substituint els caràcters '_' per '.', i generant la capçalera del fitxer FASTA amb nom de l'haplotip,
# nombre de reads i freq relativa. També es substitueix el terme Hpl per HplFw amb la funció 'gsub()'.
names(seqsf) <- paste(gsub("Hpl","HplFw",gsub('_','.',names(seqsf))),nrf,frq_lst,sep="|")
## Guarda a la taula de resultats per a cadenes fw:
# El total de reads de la mostra avaluada
rdf.fw$fw.all[i] <- sum(lstf$nr)
# El nº de reads amb baixa freqüència (menor al mínim permès)
rdf.fw$fw.lowf[i] <- sum(lstf$nr)-sum(nrf)
## Aplica el mateix procés per a les cadenes classificades reverse de la mostra avaluada
# Aplica la funció 'ReadAmplSeqs()' del paquet QSutils que permet obtenir els haplotips i les
# seves freqüències a partir d'un fitxer FASTA indicat.
lstr <- ReadAmplSeqs(trimfiles[idx.rv[i]])
# Guarda les seqüències de longitud major al mínim permès i filtra el nº de reads
# eliminant els de les seqüències que presenten longitud menor al mínim permès
lstr$nr <- lstr$nr[nchar(lstr$hseqs) >= min.seq.len]
# Aplica la funció 'GetQSData()' del paquet QSutils per llegir els amplicons del
# fitxer amb els valores d'abundància, filtrar per mínima abundància i ordenar
# segons les mutacions.
s2 <- GetQSData(trimfiles[idx.rv[i]],thr)
# Guarda les seqüències de longitud major al mínim permès
rfilt <- nchar(s2$seqs) >= min.seq.len
# Filtra les seqüències per eliminar les que presenten longitud menor al mínim permès
seqsr <- s2$seqs[rfilt]
nrr <- s2$nr[rfilt]
# Calcula la freqüència relativa de cada amplicó filtrat per abundància respecte el total
# de reads del fitxer inicial.
frq_lst <- round(nrr/sum(lstr$nr)*100,2)
# Donat que la funció del paquet QSutils retorna els noms dels haplotips diferent al desitjat, es canvien els
# noms substituint els caràcters '_' per '.', i generant la capçalera del fitxer FASTA amb nom de l'haplotip,
# nombre de reads i freq relativa. També es substitueix el terme Hpl per HplFw amb la funció 'gsub()'
names(seqsr) <- paste(gsub("Hpl","HplRv",gsub('_','.',names(seqsr))),nrr,frq_lst,sep="|")
## Guarda a la taula de resultats per a cadenes rv:
# El total de reads de la mostra avaluada
rdf.rv$rv.all[i] <- sum(lstr$nr)
# El nº de reads amb baixa freqüència (menor al mínim permès)
rdf.rv$rv.lowf[i] <- sum(lstr$nr)-sum(nrr)
# Guarda en un vector les seqs de les dues cadenes
aseqs <- c(seqsf,seqsr)
### Aliniament múltiple per muscle dels haplotips fw i rv
# Guarda el resultat de l'aliniament múltiple realitzat amb la funció 'muscle()'
ma <- DNAStringSet(muscle::muscle(aseqs,log = "./tmp/muscle.log", verbose= TRUE,quiet=TRUE))
# Genera el fitxer resultant de l'aliniament múltiple al directori MACH per guardar l'aliniament
# dels haplotips de la mostra avaluada
writeXStringSet(ma,ma.flnms[i])
# Aplica la funció per obtenir els noms dels haplotips i les seves dades
lst <- ReadAmplSeqs(ma.flnms[i])
# Guarda el vector que inclou el nº de seqüències dels haplotips
nr <- lst$nr
# Guarda les seqüències dels haplotips amb més reads del mínim permès
ma_seqs <- lst$hseqs
## Separar seq FW i RV ja aliniades
# Després de l'aliniament múltiple de totes les seqüències de la mostra, separa en 2 variables les que
# corresponen a cadena forward o reverse
ifw <- grep("^HplFw",names(ma_seqs))
irv <- grep("^HplRv",names(ma_seqs))
# Aplica la funció del paquet QSutils que permet calcular la intersecció entre els haplotips
lstI<-IntersectStrandHpls(nr[ifw],ma_seqs[ifw],nr[irv],ma_seqs[irv],thr=0)
# Per obtenir realment les seqs coincidents entre ambdues cadenes cal afegir el codi següent:
# A partir dels vectors d'abundància de les cadenes forward i reverse aliniades, indica
# en quins casos coincideixen les seqs, ja que seran els indexs on el nº de reads
# per FW i per RV sigui major a 0
fl <- lstI$pFW > 0 & lstI$pRV > 0
# Combina les dades dels dos conjunts de seqüències, eliminant aquelles que estiguin duplicades entre les cadenes
# (només apareixeran un cop)
nms <- union(ma_seqs[ifw],ma_seqs[irv])
# Filtra les seqüències dels haplotips coincidents en FW i RV
hseqs <- nms[fl]
# Calcula l'abundància (nº reads) dels haplotips en ambdues cadenes
nr <- lstI$pFW[fl]+lstI$pRV[fl]
# Ordena les abundàncies en ordre decreixent
o <- order(nr, decreasing = TRUE)
# Reordena les seqs en funció de la seva abundància
hseqs=hseqs[o]
# Reordena els valors d'abundància segons els indexs obtinguts
nr=nr[o]
# Calcula la freq relativa (nº de reads d'un haplotip entre el total dels que entraven a la intersecció d'aquella cadena)
# per aquells haplotips coincidents en ambdues cadenes, i indica el valor mínim de freq al comparar cada parella d'haplotips
# coincidents -> freq mínima d'intersecció
ov <- pmin((lstI$pFW/sum(lstI$pFW)),(lstI$pRV/sum(lstI$pRV)))
## Guarda a les taules de resultats de cadenes FW i RV:
# El nº total de reads dels haplotips que entraven a la intersecció
rdf.fw$fw.in[i] <- sum(lstI$pFW)
rdf.rv$rv.in[i] <- sum(lstI$pRV)
# El nº de reads dels haplotips coincidents entre ambdues cadenes
rdf.fw$fw.com[i] <- sum(lstI$pFW[fl])
rdf.rv$rv.com[i] <- sum(lstI$pRV[fl])
# El nº de reads dels haplotips no coincidents (únics d'una cadena)
rdf.fw$fw.unq[i] <- sum(lstI$pFW[!fl])
rdf.rv$rv.unq[i] <- sum(lstI$pRV[!fl])
## Guarda a la taula de resultats globals (que en aquest punt encara està buida):
# El total de reads FW+RV de la mostra avaluada
rdf.gbl$all[i] <- rdf.fw$fw.all[i]+rdf.rv$rv.all[i]
# Sumatori d reads dels haplotips amb baixa freqüència d'ambues cadenes
rdf.gbl$lowf[i] <- rdf.fw$fw.lowf[i] + rdf.rv$rv.lowf[i]
# Sumatori dels reads únics d'una cadena de DNA (no coincidents)
rdf.gbl$unq[i] <- rdf.fw$fw.unq[i] + rdf.rv$rv.unq[i]
# Sumatori de tots els valors mínims de freq calculats entre les parelles coincidents: intersecció entre reads
# Percentatge de superfície que solapa entre ambdues cadenes
rdf.gbl$ovrlp[i] <- round(sum(ov)*100,2)
# Percentatge de reads dels haplotips coincidents en FW i RV: % de reads comuns respecte el total
# Coincideix amb el resultat de dividir els reads que han coincidit (suma de fw.com i rv.com) entre els reads
# que han entrat en la intersecció (suma de fw.in i rv.in)
rdf.gbl$common[i] <- round((lstI$nr)/(rdf.fw$fw.in[i]+rdf.rv$rv.in[i])*100,2)
# Sumatori del nº de reads totals dels haplotips coincidents d'ambdues cadenes: reads finals
rdf.gbl$Fn.rd[i] <- lstI$nr
### Només guarda els resultats pertinents si detecta intersecció entre les cadenes
if(length(nr)!=0) {
### Gràfic de barres dels haplotips alineats
# Concatenació de l'ID del pacient i la regió avaluada en la iteració
tt <- paste(FlTbl$Pat.ID[idx.fw[i]]," - ",FlTbl$Ampl.Nm[idx.fw[i]],sep="")
# Suma el nº de seqs dels haplotips FW i RV independentment de si coincideixen o no
p <- lstI$pFW+lstI$pRV
# Substitueix el nº de reads d'aquells haplotips no coincidents per 0: Només queda el sumatori de reads
# per cada haplotip que ha coincidit en ambdues cadenes
p[ lstI$pFW==0 | lstI$pRV==0 ] <- 0
# Calcula la freq relativa del nº de reads de cada haplotip entre el total dels que han coincidit
p <- p/sum(p)
## Genera els gràfics per representar els haplotips de cada cadena i els aliniats amb les seves freq
# Divideix el nº de seqs de cada haplotip de la cadena FW entre el total dels que entraven a la intersecció
pA <- lstI$pFW/sum(lstI$pFW)
# Divideix el nº de seqs de cada haplotip de la cadena RV entre el total
pB <- lstI$pRV/sum(lstI$pRV)
# Defineix el límit de l'eix Y en funció de les freqüències calculades
ymx <- max(c(pA,pB))
# Gràfic de barres per representar els haplotips de la cadena FW i la seva freq relativa
barplot(pA,ylim=c(0,ymx)); abline(h=0)
title(main="FW strand haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
# Gràfic de barres per representar els haplotips de la cadena RV i la seva freq relativa
barplot(pB,ylim=c(0,ymx)); abline(h=0)
title(main="RV strand haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
# Gràfic de barres per representar els haplotips coincidents en ambdues cadenes i la seva freq
barplot(p); abline(h=0)
title(main="Intersected haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
### Aplica la funció del mateix paquet per guardar els haplotips coincidents en ambdues cadenes
# amb les seves freqüències en els fitxers MACHpl02
s<-SaveHaplotypes(out.flnms[i],as.character(hseqs),nr=nr,max.difs=max.difs)
}
### Si no es detecta superposició entre els haplotips d'ambdues cadenes, salta a la seqüent iteració
else if(length(nr)==0) {
cat("\n--------------------------------------------------\n")}
} # Fi del loop sobre totes les mostres (2 per pacient)
dev.off()
# Genera el fitxer .txt de resultats de les interseccions
sink(file=file.path(repDir,"MA.Intersects-SummRprt.txt"))
cat("\n FW + RV HAPLOTYPES INTERSECTIONS")
cat("\n======================================\n")
cat("\nCutting FW and RV at ",thr,"% followed by haplotypes intersection.\n",
sep="")
cat("\nFrequencies as sum of reads FW+RV.\n\n")
## Resultats per total de reads
cat("\nSUMMARY RESULTS BY READS NUMBER")
cat("\n===============================\n\n")
# A partir de la taula FlTbl (provinent de l'arxiu RData del pas anterior del pipeline), guarda
# les entrades corresponents a la cadena forward
fl.fw <- FlTbl$Str=="fw"
# Genera un data frame amb els ID dels pacients, la regió amplificada, i els resultats obtinguts
# al llarg de la intersecció d'haplotips per la cadena FW
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.fw,stringsAsFactors=FALSE)
print(frdf)
cat("\n")
# Genera el mateix data frame d'abans pels haplotips de la cadena RV per incloure-la al fitxer
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.rv,stringsAsFactors=FALSE)
print(frdf)
cat("\n")
# Genera un altre data frame per incloure els resultats globals d'ambdues cadenes (FW+RV)
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.gbl,stringsAsFactors=FALSE)
print(frdf)
## Mostra els resultats d'aplicar el sumatori sobre les columnes de resultats
# de cadascun dels 3 data frames anteriors
cat("\n\nTotal counts:\n\n")
tots.fw <- apply(rdf.fw,2,sum)
print(tots.fw)
cat("\n")
tots.rv <- apply(rdf.rv,2,sum)
print(tots.rv)
cat("\n")
tots.gbl <- apply(rdf.gbl[,c(1:3,6)],2,sum)
print(tots.gbl)
cat("\n")
## Aplica els percentatges sobre els resultats dels sumatoris anterior per FW, RV y global
cat("\nPercentage:\n\n")
ptts <- round(tots.fw/tots.fw[1]*100,2)
print(ptts)
cat("\n")
ptts <- round(tots.rv/tots.rv[1]*100,2)
print(ptts)
cat("\n")
ptts <- round(tots.gbl/tots.gbl[1]*100,2)
print(ptts)
cat("\n")
sink()
## Genera el pdf que es guardarà a la carpeta reports on s'inclouran els gràfics
# de resultats globals de la intersecció d'haplotips
pdf.flnm2 <- "IntersectBarplots.pdf"
pdf(file.path(repDir,pdf.flnm2),paper="a4r",width=10,height=6)
par(mar=c(7,4,4,2)+.1)
## Genera la paleta de colors
cls <- brewer.pal(8,"Dark2")
pal <- cls[1:2]
# Guarda els noms de les mostres a partir de la concatenació de l'ID del pacient i
# la regió avaluada (amb el condicional FW per no repetir els noms dos cops)
bnms <- paste(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw])
# Genera una taula amb les columnes corresponents als valors de reads comuns en FW i RV
vals <- cbind(rdf.fw$fw.com,rdf.rv$rv.com)
# Defineix el límit màxim de l'eix Y a partir del sumatori de reads comuns (aplicat sobre files)
ymx <- max(rowSums(vals))*1.2
# Gràfic de barres on es representen els reads comuns en ambdues cadenes en funció de la mostra
# names.arg correspon a un vector on s'indiquen els noms a representar sota les barres del gràfic
barplot(t(vals),col=pal,ylim=c(0,ymx),names.arg=bnms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("fw","rv"),cex=0.8)
title(main="Intersected reads")
## Genera una altra taula on s'inclouen les columnes de reads de baixa freqüència i únics d'una cadena
# per als haplotips FW i RV, fent el sumatori per files (per mostra)
vals <- cbind(rowSums(rdf.fw[,c("fw.lowf","fw.unq")]),
rowSums(rdf.rv[,c("rv.lowf","rv.unq")]))
# Defineix el límit màxim de l'eix Y a partir del sumatori aplicat sobre files de la taula anterior
ymx <- max(rowSums(vals))*1.2
# Gràfic de barres on es representen els reads filtrats totals (per abundància o per no trobar intersecció)
# en funció de la mostra
barplot(t(vals),col=pal,ylim=c(0,ymx),names.arg=bnms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("fw","rv"),cex=0.8)
title(main="Filtered out reads")
## Defineix la nova paleta de colors
pal <- cls[c(1,3,2,4)]
# Guarda en una matriu les columnes amb els valors de reads finals, únics de cadena i de baixa freq
# (total per ambdues cadenes)
mbp <- data.matrix(frdf[,c(8,5,4)])
# Guarda els noms de les mostres a partir de la concatenació de l'ID del pacient i
# la regió avaluada separats per punt (aquest cop els agafa de la taula de resultats)
mbp.nms <- paste(frdf$Pat.ID,frdf$Ampl.Nm,sep=".")
omar <- par(mar=c(6,3.5,3,2)) # Defineix els marges de pàgina
# Gràfic que representa, per cada mostra, una barra indicant en diferents colors els reads
# comuns (els finals després de la intersecció), els únics de cadena i els de baixa freq
barplot(t(mbp),col=pal,ylim=c(0,max(rowSums(mbp))*1.2),names.arg=mbp.nms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("Common","Unique","Low freq."),
cex=0.8)
## Genera un altre gràfic de barres per representar els percentatges de reads
# comuns per cada mostra (rendiment de la intersecció)
par(mar=omar)
barplot(frdf$common,col="lavender",ylim=c(0,100),names.arg=mbp.nms,
las=2,cex.names=0.6,cex.axis=0.8)
title(main="FW + RV intersection yield",line=1)
# Fa coerció de la columna reads per tal de poder fer el sumatori més endavant
FlTbl$Reads <- as.numeric(FlTbl$Reads)
## Guarda els noms concatenants l'ID dels pacients i la regió amplificada separats per punt
all.nms <- paste(FlTbl$Pat.ID,FlTbl$Ampl.Nm,sep=".")
# Realitza el sumatori dels reads FW i RV de cadascuna de les mostres
trds <- tapply(FlTbl$Reads,all.nms,sum)
# Actualitza el vector amb els noms per eliminar els duplicats, ja que al fer el sumatori
# pasa a tenir una única entrada per mostra (i no 2 quan tenia FW i RV)
all.nms <- names(trds)
# Guarda la columna de reads finals per mostra (després de la intersecció)
frds <- frdf$Fn.rd
# Els noms de la columna de reads finals seran de nou la concatenació de ID de pacient i regió
names(frds) <- paste(frdf$Pat.ID,frdf$Ampl.Nm,sep=".")
# Genera una matriu buida amb tantes files com noms de mostres s'han generat, i 2 columnes
ytbl <- matrix(0,nrow=length(all.nms),ncol=2)
# Els noms de les files de la matriu correspondran als noms de les mostres
rownames(ytbl) <- all.nms
# Els noms de les columnes corresponen als reads totals i els que han quedat
# després de la intersecció
colnames(ytbl) <- c("All","Passed")
# Guarda els valors de la columna de reads totals
ytbl[,1] <- trds
# Assigna a cada mostra els valors de la columna de reads finals (que han passat)
ytbl[names(frds),2] <- frds
omar <- par(mar=c(6,3.5,3,2)) # Marges de pàgina
# Gràfic de barres que representa, per cada mostra, una barra indicant els reads totals
# i una altra els reads que han passat
barplot(t(ytbl),beside=TRUE,col=cls[1:2],ylim=c(0,max(ytbl)*1.2),las=2,
names.arg=rownames(ytbl),cex.names=0.6,cex.axis=0.8)
abline(h=0)
legend("top",horiz=TRUE,fill=cls,legend=colnames(ytbl),cex=0.8)
## Últim gràfic per representar el rendiment de la intersecció, dividint els
# reads que han passat respecte els totals (els inicials del pas)
par(mar=omar)
barplot(as.vector(ytbl[,2]/ytbl[,1]*100),col="lavender",ylim=c(0,100),
names.arg=rownames(ytbl),las=2,cex.names=0.6,cex.axis=0.8)
abline(h=0)
title(main="Global yield",line=1)
dev.off()
# Retorna el data frame amb els resultats globals
return(frdf)
}
aliafdz/QApckg documentation built on June 2, 2022, 10:29 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ConsHaplotypes
Documented in ConsHaplotypes
#' @title Generate consensus haplotypes
#'
#' @description Computes the intersection of forward and reverse strand haplotypes and generates some report files.
#'
#' @details This function is designed to be used after the execution of \code{\link{demultiplexPrimer}} function
#' from the same package. After the generation of FASTA files containing forward and reverse strand reads
#' for the evaluated samples, \code{\link{ConsHaplotypes}} executes multiple alignment with \code{\link{muscle}} and
#' returns the consensus haplotypes using \code{\link{IntersectStrandHpls}}, that will be saved using
#' the helper function \code{\link{SaveHaplotypes}}.
#'
#'
#' @param trimfiles Vector including the paths of demultiplexed files by specific primer, with fna extension.
#' @param pm.res The list returned by \code{\link{demultiplexPrimer}}, including \code{fileTable}
#' and \code{poolTable} data frames.
#' @param thr Threshold to filter haplotypes at minimum abundance before multiple alignment.
#' @param min.seq.len Threshold to filter haplotypes at minimum length before intersection.
#' @param max.difs Maximum number of mismatches allowed in resulting consensus haplotypes with
#' respect to the dominant one.
#'
#' @return The function returns a \code{\link{data.frame}} object containing the intersection results
#' for each combination of patient and amplicon region, including the initial number of reads, filtered out reads
#' (for being below a given frequency threshold or unique to a single strand), overlapping frequency
#' between both strands and the common reads (in percentage and nº of reads).
#'
#' @return After execution, two FASTA files for each combination of sample and pool will be saved in a newly
#' generated MACH folder; the first includes multiple alignment between forward and reverse strand haplotypes,
#' and the second includes the forward and reverse strands intersected.
#' Additionaly, some report files will be generated in the reports folder:
#' \enumerate{
#' \item{\code{MA.Intersects-SummRprt.txt}: Includes the sumary results by reads number after abundance filter and
#' strand intersection.}
#' \item{\code{MA.Intersects.plots.pdf}: Includes different barplots for each sample representing the frequency of
#' forward, reverse and intersected strand haplotypes.}
#' \item{\code{IntersectBarplots.pdf}: Includes different barplots for all combinations of patient and pool,
#' representing the number of intersected and filtered out reads, the intersection yield and global yield.}}
#'
#' @note A new file named \code{muscle.log} containing \code{\link{muscle}} options will be generated and saved in a folder named "tmp".
#'
#' @import Biostrings
#' @import ape
#' @import QSutils
#' @import parallel
#' @importFrom muscle muscle
#' @importFrom RColorBrewer brewer.pal
#' @importFrom methods as
#' @importFrom methods is
#' @seealso \code{\link{muscle}}, \code{\link{IntersectStrandHpls}}, \code{\link{demultiplexPrimer}}, \code{\link{SaveHaplotypes}}
#' @examples
#' splitDir <- "./splits"
#' # Save the file names with complete path
#' splitfiles <- list.files(splitDir,recursive=TRUE,full.names=TRUE,include.dirs=TRUE)
#' # Get data
#' samples <- read.table("./data/samples.csv", sep="\t", header=T,
#' colClasses="character",stringsAsFactors=F)
#' mids <- read.table("./data/mids.csv", sep="\t", header=T,
#' stringsAsFactors=F)
#' # Apply previous function from QA analysis
#' pm.res <- demultiplexPrimer(splitfiles,samples,primers)
#' # Save the files generated by previous function
#' trimDir <- "./trim"
#' trimfiles <- list.files(trimDir,recursive=TRUE,full.names=TRUE,include.dirs=TRUE)
#' # Define necessary parameters
#' min.seq.len <- 150
#' thr <- 0.2
#' int.res <- ConsHaplotypes(trimfiles, pm.res, thr, min.seq.len)
#'
#' @export
#' @author Alicia Aranda
ConsHaplotypes <- function(trimfiles,pm.res,thr=0.2, min.seq.len=150,max.difs=250){
# Si la ruta on es troben els fitxers de la carpeta trim no està ben especificada, intenta buscar la
# ruta correcta a partir del directori de treball
# Si tot i així no troba fitxers, atura l'execució i mostra un missatge d'error
if(length(trimfiles)==0) {
trimfiles <- list.files(paste(getwd(),"/trim",sep=''))
if(length(trimfiles)==0) {
stop("Please indicate correct path for demultiplexed reads by primer.\n")
}
}
# Si la taula pm.res no s'ha inclòs o no existeix, retorna un missatge d'error
if(missing(pm.res)|!exists('pm.res')|length(pm.res)==0) {
stop("The list obtained by demultiplexPrimer function is needed.\n")
}
if(!is.list(pm.res)){
stop("pm.res argument must be a list.\n")
}
# Si cap dels fitxers indicats a la carpeta trim no existeix, atura l'execució i
# mostra un missatge d'error
if(any(!file.exists(trimfiles))){
stop(paste(trimfiles[!file.exists(trimfiles)],"does not exist.\n"))
}
# Si no existeixen les carpetes on es guarden els fitxers resultants de la funció,
# es generen automàticament a la carpeta de treball
if(!dir.exists("./MACH")) {
dir.create("./MACH")}
mach.Dir <- "./MACH"
if(!dir.exists("./reports")) {
dir.create("./reports")}
repDir <- "./reports"
if(!dir.exists("./tmp")) {
dir.create("./tmp")}
# En lloc de carregar el RData es carreguen les taules fileTable i poolTable de la variable de la funció anterior.
FlTbl<- pm.res$fileTable
PoolTbl<- pm.res$poolTable
# Retorna els indexs de la taula FlTbl derivats d'ordenar les mostres en funció de l'ID del
# pacient, de l'amplicó (regió) avaluat i la cadena forward o reverse
o <- order(FlTbl$Pat.ID,FlTbl$Ampl.Nm,FlTbl$Str)
# Reordena les entrades de la taula en funció dels indexs anteriors, de manera que s'agrupen
# les entrades corresponents al mateix pacient i en segon ordre a la mateixa regió avaluada
FlTbl <- FlTbl[o,]
### Indica quins dels fitxers FASTA inclosos al directori trim corresponen a
# cadascuna de les cadenes
idx.fw<-grep("PrFW",trimfiles)
idx.rv<-grep("PrRV",trimfiles)
# Guarda els noms dels fitxers resultants d'aquest script, que correspondran a la concatenació
# del terme MACHpl02, l'ID del pacient i el nom de la regió avaluada, indicats de manera que només
# s'obtingui un fitxer .fna per a cada regió del mateix pacient (per això el condicional fw)
out.flnms <- paste("MACHpl02",FlTbl$Pat.ID[idx.fw],
FlTbl$Ampl.Nm[idx.fw],"fna",sep=".")
# Guarda la ruta on es desaran els fitxers .fna, a la carpeta MACH
out.flnms <- file.path(mach.Dir,out.flnms)
# Guarda els noms d'uns altres fitxers resultants, que correspondran a la concatenació
# del terme MAfwrv, l'ID del pacient i el nom de la regió avaluada, indicats de manera que només
# s'obtingui un fitxer .fna per a cada regió del mateix pacient (per això el condicional fw)
ma.flnms <- paste("MAfwrv",FlTbl$Pat.ID[idx.fw],
FlTbl$Ampl.Nm[idx.fw],"fna",sep=".")
# Guarda la ruta on es desaran els fitxers .fna, a la carpeta MACH
ma.flnms <- file.path(mach.Dir,ma.flnms)
# Guarda la meitat de la longitud de fitxers presents a la carpeta trim del fitxer previ
n <- length(trimfiles)/2
# Genera un data frame inicialment buit, amb n files (2 entrades per pacient) i 5 columnes
# que inclouran els resultats obtinguts en pasos posteriors per a les cadenes fw
rdf.fw <- data.frame(fw.all=integer(n),fw.lowf=integer(n),
fw.in=integer(n),fw.unq=integer(n),fw.com=integer(n))
# Genera un altre data frame també buit, amb n files i 5 columnes per incloure resultats
# de les cadenes rv
rdf.rv <- data.frame(rv.all=integer(n),rv.lowf=integer(n),
rv.in=integer(n),rv.unq=integer(n),rv.com=integer(n))
# Genera un altre data frame buit, amb n files i 6 columnes, per incloure els resultats globals
rdf.gbl <- data.frame(all=integer(n),lowf=integer(n),unq=integer(n),
ovrlp=numeric(n),common=numeric(n),Fn.rd=integer(n))
# Genera el primer fitxer pdf on s'inclouran diverses representacions de resultats generats en el bucle for
# Per cada regió amplificada de cada pacient, es representen els haplotips de cada cadena i els que han
# coincidit entre les dues, amb les corresponents freqüències
pdf(file.path(repDir,"MA.Intersects.plots.pdf"),paper="a4",
width=7,height=11)
par(mfrow=c(3,1))
# Bucle per iterar sobre totes les mostres (2 per pacient)
foreach(i=c(1:n)) %do% {
# Aplica la funció 'ReadAmplSeqs()' del paquet QSutils que permet obtenir els haplotips i les
# seves freqüències a partir d'un fitxer FASTA indicat.
lstf <- ReadAmplSeqs(trimfiles[idx.fw[i]])
# Guarda les seqüències de longitud major al mínim permès i filtra el nº de reads
# eliminant els de les seqüències que presenten longitud menor al mínim permès
lstf$nr <- lstf$nr[nchar(lstf$hseqs) >= min.seq.len]
# Aplica la funció 'GetQSData()' del paquet QSutils per llegir els amplicons del
# fitxer amb els valores d'abundància, filtrar per mínima abundància i ordenar
# segons les mutacions.
s1 <- GetQSData(trimfiles[idx.fw[i]],thr)
# Guarda les seqüències de longitud major al mínim permès
ffilt <- nchar(s1$seqs) >= min.seq.len
# Filtra les seqüències per eliminar les que presenten longitud menor al mínim permès
seqsf <- s1$seqs[ffilt]
nrf <- s1$nr[ffilt]
# Calcula la freqüència relativa de cada amplicó filtrat per abundància respecte el total
# de reads del fitxer inicial.
frq_lst <- round(nrf/sum(lstf$nr)*100,2)
# Donat que la funció del paquet QSutils retorna els noms dels haplotips diferent al desitjat, es canvien els
# noms substituint els caràcters '_' per '.', i generant la capçalera del fitxer FASTA amb nom de l'haplotip,
# nombre de reads i freq relativa. També es substitueix el terme Hpl per HplFw amb la funció 'gsub()'.
names(seqsf) <- paste(gsub("Hpl","HplFw",gsub('_','.',names(seqsf))),nrf,frq_lst,sep="|")
## Guarda a la taula de resultats per a cadenes fw:
# El total de reads de la mostra avaluada
rdf.fw$fw.all[i] <- sum(lstf$nr)
# El nº de reads amb baixa freqüència (menor al mínim permès)
rdf.fw$fw.lowf[i] <- sum(lstf$nr)-sum(nrf)
## Aplica el mateix procés per a les cadenes classificades reverse de la mostra avaluada
# Aplica la funció 'ReadAmplSeqs()' del paquet QSutils que permet obtenir els haplotips i les
# seves freqüències a partir d'un fitxer FASTA indicat.
lstr <- ReadAmplSeqs(trimfiles[idx.rv[i]])
# Guarda les seqüències de longitud major al mínim permès i filtra el nº de reads
# eliminant els de les seqüències que presenten longitud menor al mínim permès
lstr$nr <- lstr$nr[nchar(lstr$hseqs) >= min.seq.len]
# Aplica la funció 'GetQSData()' del paquet QSutils per llegir els amplicons del
# fitxer amb els valores d'abundància, filtrar per mínima abundància i ordenar
# segons les mutacions.
s2 <- GetQSData(trimfiles[idx.rv[i]],thr)
# Guarda les seqüències de longitud major al mínim permès
rfilt <- nchar(s2$seqs) >= min.seq.len
# Filtra les seqüències per eliminar les que presenten longitud menor al mínim permès
seqsr <- s2$seqs[rfilt]
nrr <- s2$nr[rfilt]
# Calcula la freqüència relativa de cada amplicó filtrat per abundància respecte el total
# de reads del fitxer inicial.
frq_lst <- round(nrr/sum(lstr$nr)*100,2)
# Donat que la funció del paquet QSutils retorna els noms dels haplotips diferent al desitjat, es canvien els
# noms substituint els caràcters '_' per '.', i generant la capçalera del fitxer FASTA amb nom de l'haplotip,
# nombre de reads i freq relativa. També es substitueix el terme Hpl per HplFw amb la funció 'gsub()'
names(seqsr) <- paste(gsub("Hpl","HplRv",gsub('_','.',names(seqsr))),nrr,frq_lst,sep="|")
## Guarda a la taula de resultats per a cadenes rv:
# El total de reads de la mostra avaluada
rdf.rv$rv.all[i] <- sum(lstr$nr)
# El nº de reads amb baixa freqüència (menor al mínim permès)
rdf.rv$rv.lowf[i] <- sum(lstr$nr)-sum(nrr)
# Guarda en un vector les seqs de les dues cadenes
aseqs <- c(seqsf,seqsr)
### Aliniament múltiple per muscle dels haplotips fw i rv
# Guarda el resultat de l'aliniament múltiple realitzat amb la funció 'muscle()'
ma <- DNAStringSet(muscle::muscle(aseqs,log = "./tmp/muscle.log", verbose= TRUE,quiet=TRUE))
# Genera el fitxer resultant de l'aliniament múltiple al directori MACH per guardar l'aliniament
# dels haplotips de la mostra avaluada
writeXStringSet(ma,ma.flnms[i])
# Aplica la funció per obtenir els noms dels haplotips i les seves dades
lst <- ReadAmplSeqs(ma.flnms[i])
# Guarda el vector que inclou el nº de seqüències dels haplotips
nr <- lst$nr
# Guarda les seqüències dels haplotips amb més reads del mínim permès
ma_seqs <- lst$hseqs
## Separar seq FW i RV ja aliniades
# Després de l'aliniament múltiple de totes les seqüències de la mostra, separa en 2 variables les que
# corresponen a cadena forward o reverse
ifw <- grep("^HplFw",names(ma_seqs))
irv <- grep("^HplRv",names(ma_seqs))
# Aplica la funció del paquet QSutils que permet calcular la intersecció entre els haplotips
lstI<-IntersectStrandHpls(nr[ifw],ma_seqs[ifw],nr[irv],ma_seqs[irv],thr=0)
# Per obtenir realment les seqs coincidents entre ambdues cadenes cal afegir el codi següent:
# A partir dels vectors d'abundància de les cadenes forward i reverse aliniades, indica
# en quins casos coincideixen les seqs, ja que seran els indexs on el nº de reads
# per FW i per RV sigui major a 0
fl <- lstI$pFW > 0 & lstI$pRV > 0
# Combina les dades dels dos conjunts de seqüències, eliminant aquelles que estiguin duplicades entre les cadenes
# (només apareixeran un cop)
nms <- union(ma_seqs[ifw],ma_seqs[irv])
# Filtra les seqüències dels haplotips coincidents en FW i RV
hseqs <- nms[fl]
# Calcula l'abundància (nº reads) dels haplotips en ambdues cadenes
nr <- lstI$pFW[fl]+lstI$pRV[fl]
# Ordena les abundàncies en ordre decreixent
o <- order(nr, decreasing = TRUE)
# Reordena les seqs en funció de la seva abundància
hseqs=hseqs[o]
# Reordena els valors d'abundància segons els indexs obtinguts
nr=nr[o]
# Calcula la freq relativa (nº de reads d'un haplotip entre el total dels que entraven a la intersecció d'aquella cadena)
# per aquells haplotips coincidents en ambdues cadenes, i indica el valor mínim de freq al comparar cada parella d'haplotips
# coincidents -> freq mínima d'intersecció
ov <- pmin((lstI$pFW/sum(lstI$pFW)),(lstI$pRV/sum(lstI$pRV)))
## Guarda a les taules de resultats de cadenes FW i RV:
# El nº total de reads dels haplotips que entraven a la intersecció
rdf.fw$fw.in[i] <- sum(lstI$pFW)
rdf.rv$rv.in[i] <- sum(lstI$pRV)
# El nº de reads dels haplotips coincidents entre ambdues cadenes
rdf.fw$fw.com[i] <- sum(lstI$pFW[fl])
rdf.rv$rv.com[i] <- sum(lstI$pRV[fl])
# El nº de reads dels haplotips no coincidents (únics d'una cadena)
rdf.fw$fw.unq[i] <- sum(lstI$pFW[!fl])
rdf.rv$rv.unq[i] <- sum(lstI$pRV[!fl])
## Guarda a la taula de resultats globals (que en aquest punt encara està buida):
# El total de reads FW+RV de la mostra avaluada
rdf.gbl$all[i] <- rdf.fw$fw.all[i]+rdf.rv$rv.all[i]
# Sumatori d reads dels haplotips amb baixa freqüència d'ambues cadenes
rdf.gbl$lowf[i] <- rdf.fw$fw.lowf[i] + rdf.rv$rv.lowf[i]
# Sumatori dels reads únics d'una cadena de DNA (no coincidents)
rdf.gbl$unq[i] <- rdf.fw$fw.unq[i] + rdf.rv$rv.unq[i]
# Sumatori de tots els valors mínims de freq calculats entre les parelles coincidents: intersecció entre reads
# Percentatge de superfície que solapa entre ambdues cadenes
rdf.gbl$ovrlp[i] <- round(sum(ov)*100,2)
# Percentatge de reads dels haplotips coincidents en FW i RV: % de reads comuns respecte el total
# Coincideix amb el resultat de dividir els reads que han coincidit (suma de fw.com i rv.com) entre els reads
# que han entrat en la intersecció (suma de fw.in i rv.in)
rdf.gbl$common[i] <- round((lstI$nr)/(rdf.fw$fw.in[i]+rdf.rv$rv.in[i])*100,2)
# Sumatori del nº de reads totals dels haplotips coincidents d'ambdues cadenes: reads finals
rdf.gbl$Fn.rd[i] <- lstI$nr
### Només guarda els resultats pertinents si detecta intersecció entre les cadenes
if(length(nr)!=0) {
### Gràfic de barres dels haplotips alineats
# Concatenació de l'ID del pacient i la regió avaluada en la iteració
tt <- paste(FlTbl$Pat.ID[idx.fw[i]]," - ",FlTbl$Ampl.Nm[idx.fw[i]],sep="")
# Suma el nº de seqs dels haplotips FW i RV independentment de si coincideixen o no
p <- lstI$pFW+lstI$pRV
# Substitueix el nº de reads d'aquells haplotips no coincidents per 0: Només queda el sumatori de reads
# per cada haplotip que ha coincidit en ambdues cadenes
p[ lstI$pFW==0 | lstI$pRV==0 ] <- 0
# Calcula la freq relativa del nº de reads de cada haplotip entre el total dels que han coincidit
p <- p/sum(p)
## Genera els gràfics per representar els haplotips de cada cadena i els aliniats amb les seves freq
# Divideix el nº de seqs de cada haplotip de la cadena FW entre el total dels que entraven a la intersecció
pA <- lstI$pFW/sum(lstI$pFW)
# Divideix el nº de seqs de cada haplotip de la cadena RV entre el total
pB <- lstI$pRV/sum(lstI$pRV)
# Defineix el límit de l'eix Y en funció de les freqüències calculades
ymx <- max(c(pA,pB))
# Gràfic de barres per representar els haplotips de la cadena FW i la seva freq relativa
barplot(pA,ylim=c(0,ymx)); abline(h=0)
title(main="FW strand haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
# Gràfic de barres per representar els haplotips de la cadena RV i la seva freq relativa
barplot(pB,ylim=c(0,ymx)); abline(h=0)
title(main="RV strand haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
# Gràfic de barres per representar els haplotips coincidents en ambdues cadenes i la seva freq
barplot(p); abline(h=0)
title(main="Intersected haplotypes barplot",line=0.5,cex.main=1)
title(sub=tt,line=0.5,cex.main=1)
### Aplica la funció del mateix paquet per guardar els haplotips coincidents en ambdues cadenes
# amb les seves freqüències en els fitxers MACHpl02
s<-SaveHaplotypes(out.flnms[i],as.character(hseqs),nr=nr,max.difs=max.difs)
}
### Si no es detecta superposició entre els haplotips d'ambdues cadenes, salta a la seqüent iteració
else if(length(nr)==0) {
cat("\n--------------------------------------------------\n")}
} # Fi del loop sobre totes les mostres (2 per pacient)
dev.off()
# Genera el fitxer .txt de resultats de les interseccions
sink(file=file.path(repDir,"MA.Intersects-SummRprt.txt"))
cat("\n FW + RV HAPLOTYPES INTERSECTIONS")
cat("\n======================================\n")
cat("\nCutting FW and RV at ",thr,"% followed by haplotypes intersection.\n",
sep="")
cat("\nFrequencies as sum of reads FW+RV.\n\n")
## Resultats per total de reads
cat("\nSUMMARY RESULTS BY READS NUMBER")
cat("\n===============================\n\n")
# A partir de la taula FlTbl (provinent de l'arxiu RData del pas anterior del pipeline), guarda
# les entrades corresponents a la cadena forward
fl.fw <- FlTbl$Str=="fw"
# Genera un data frame amb els ID dels pacients, la regió amplificada, i els resultats obtinguts
# al llarg de la intersecció d'haplotips per la cadena FW
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.fw,stringsAsFactors=FALSE)
print(frdf)
cat("\n")
# Genera el mateix data frame d'abans pels haplotips de la cadena RV per incloure-la al fitxer
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.rv,stringsAsFactors=FALSE)
print(frdf)
cat("\n")
# Genera un altre data frame per incloure els resultats globals d'ambdues cadenes (FW+RV)
frdf <- data.frame(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw],
rdf.gbl,stringsAsFactors=FALSE)
print(frdf)
## Mostra els resultats d'aplicar el sumatori sobre les columnes de resultats
# de cadascun dels 3 data frames anteriors
cat("\n\nTotal counts:\n\n")
tots.fw <- apply(rdf.fw,2,sum)
print(tots.fw)
cat("\n")
tots.rv <- apply(rdf.rv,2,sum)
print(tots.rv)
cat("\n")
tots.gbl <- apply(rdf.gbl[,c(1:3,6)],2,sum)
print(tots.gbl)
cat("\n")
## Aplica els percentatges sobre els resultats dels sumatoris anterior per FW, RV y global
cat("\nPercentage:\n\n")
ptts <- round(tots.fw/tots.fw[1]*100,2)
print(ptts)
cat("\n")
ptts <- round(tots.rv/tots.rv[1]*100,2)
print(ptts)
cat("\n")
ptts <- round(tots.gbl/tots.gbl[1]*100,2)
print(ptts)
cat("\n")
sink()
## Genera el pdf que es guardarà a la carpeta reports on s'inclouran els gràfics
# de resultats globals de la intersecció d'haplotips
pdf.flnm2 <- "IntersectBarplots.pdf"
pdf(file.path(repDir,pdf.flnm2),paper="a4r",width=10,height=6)
par(mar=c(7,4,4,2)+.1)
## Genera la paleta de colors
cls <- brewer.pal(8,"Dark2")
pal <- cls[1:2]
# Guarda els noms de les mostres a partir de la concatenació de l'ID del pacient i
# la regió avaluada (amb el condicional FW per no repetir els noms dos cops)
bnms <- paste(Pat.ID=FlTbl$Pat.ID[fl.fw],Ampl.Nm=FlTbl$Ampl.Nm[fl.fw])
# Genera una taula amb les columnes corresponents als valors de reads comuns en FW i RV
vals <- cbind(rdf.fw$fw.com,rdf.rv$rv.com)
# Defineix el límit màxim de l'eix Y a partir del sumatori de reads comuns (aplicat sobre files)
ymx <- max(rowSums(vals))*1.2
# Gràfic de barres on es representen els reads comuns en ambdues cadenes en funció de la mostra
# names.arg correspon a un vector on s'indiquen els noms a representar sota les barres del gràfic
barplot(t(vals),col=pal,ylim=c(0,ymx),names.arg=bnms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("fw","rv"),cex=0.8)
title(main="Intersected reads")
## Genera una altra taula on s'inclouen les columnes de reads de baixa freqüència i únics d'una cadena
# per als haplotips FW i RV, fent el sumatori per files (per mostra)
vals <- cbind(rowSums(rdf.fw[,c("fw.lowf","fw.unq")]),
rowSums(rdf.rv[,c("rv.lowf","rv.unq")]))
# Defineix el límit màxim de l'eix Y a partir del sumatori aplicat sobre files de la taula anterior
ymx <- max(rowSums(vals))*1.2
# Gràfic de barres on es representen els reads filtrats totals (per abundància o per no trobar intersecció)
# en funció de la mostra
barplot(t(vals),col=pal,ylim=c(0,ymx),names.arg=bnms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("fw","rv"),cex=0.8)
title(main="Filtered out reads")
## Defineix la nova paleta de colors
pal <- cls[c(1,3,2,4)]
# Guarda en una matriu les columnes amb els valors de reads finals, únics de cadena i de baixa freq
# (total per ambdues cadenes)
mbp <- data.matrix(frdf[,c(8,5,4)])
# Guarda els noms de les mostres a partir de la concatenació de l'ID del pacient i
# la regió avaluada separats per punt (aquest cop els agafa de la taula de resultats)
mbp.nms <- paste(frdf$Pat.ID,frdf$Ampl.Nm,sep=".")
omar <- par(mar=c(6,3.5,3,2)) # Defineix els marges de pàgina
# Gràfic que representa, per cada mostra, una barra indicant en diferents colors els reads
# comuns (els finals després de la intersecció), els únics de cadena i els de baixa freq
barplot(t(mbp),col=pal,ylim=c(0,max(rowSums(mbp))*1.2),names.arg=mbp.nms,
las=2,cex.names=0.6,cex.axis=0.8)
legend("top",horiz=TRUE,fill=pal,legend=c("Common","Unique","Low freq."),
cex=0.8)
## Genera un altre gràfic de barres per representar els percentatges de reads
# comuns per cada mostra (rendiment de la intersecció)
par(mar=omar)
barplot(frdf$common,col="lavender",ylim=c(0,100),names.arg=mbp.nms,
las=2,cex.names=0.6,cex.axis=0.8)
title(main="FW + RV intersection yield",line=1)
# Fa coerció de la columna reads per tal de poder fer el sumatori més endavant
FlTbl$Reads <- as.numeric(FlTbl$Reads)
## Guarda els noms concatenants l'ID dels pacients i la regió amplificada separats per punt
all.nms <- paste(FlTbl$Pat.ID,FlTbl$Ampl.Nm,sep=".")
# Realitza el sumatori dels reads FW i RV de cadascuna de les mostres
trds <- tapply(FlTbl$Reads,all.nms,sum)
# Actualitza el vector amb els noms per eliminar els duplicats, ja que al fer el sumatori
# pasa a tenir una única entrada per mostra (i no 2 quan tenia FW i RV)
all.nms <- names(trds)
# Guarda la columna de reads finals per mostra (després de la intersecció)
frds <- frdf$Fn.rd
# Els noms de la columna de reads finals seran de nou la concatenació de ID de pacient i regió
names(frds) <- paste(frdf$Pat.ID,frdf$Ampl.Nm,sep=".")
# Genera una matriu buida amb tantes files com noms de mostres s'han generat, i 2 columnes
ytbl <- matrix(0,nrow=length(all.nms),ncol=2)
# Els noms de les files de la matriu correspondran als noms de les mostres
rownames(ytbl) <- all.nms
# Els noms de les columnes corresponen als reads totals i els que han quedat
# després de la intersecció
colnames(ytbl) <- c("All","Passed")
# Guarda els valors de la columna de reads totals
ytbl[,1] <- trds
# Assigna a cada mostra els valors de la columna de reads finals (que han passat)
ytbl[names(frds),2] <- frds
omar <- par(mar=c(6,3.5,3,2)) # Marges de pàgina
# Gràfic de barres que representa, per cada mostra, una barra indicant els reads totals
# i una altra els reads que han passat
barplot(t(ytbl),beside=TRUE,col=cls[1:2],ylim=c(0,max(ytbl)*1.2),las=2,
names.arg=rownames(ytbl),cex.names=0.6,cex.axis=0.8)
abline(h=0)
legend("top",horiz=TRUE,fill=cls,legend=colnames(ytbl),cex=0.8)
## Últim gràfic per representar el rendiment de la intersecció, dividint els
# reads que han passat respecte els totals (els inicials del pas)
par(mar=omar)
barplot(as.vector(ytbl[,2]/ytbl[,1]*100),col="lavender",ylim=c(0,100),
names.arg=rownames(ytbl),las=2,cex.names=0.6,cex.axis=0.8)
abline(h=0)
title(main="Global yield",line=1)
dev.off()
# Retorna el data frame amb els resultats globals
return(frdf)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.