demuxStatsFromXML.R

#!/usr/bin/env Rscript

# R version : 4.0.4
## module load system/R-4.0.4_gcc-9.3.0

# demuxStatsFromXML.R
# Lecture d'un fichier XML pour extraction et mise en forme des statistiques de démultiplexage (orienté 10X pour le moment)
# Par échantillon, ce script récupère tous les index associés, le nombre de reads trouvés, dont le nombre de barcodes lus parfaitement et le nombre de barcode lus avec un mismatch.
# Ce sctipt récupère aussi les index très souvent retrouvés mais non associé à un echantillon
# Le pourcentage du nombre de fragments par échantillon sur le nombre total est calculé

## --------------------
#       PACKAGES
## --------------------
library('xml2')
library('stringr')
library('optparse')

## --------------------
#       FUNCTIONS
## --------------------
concat_df = function(df1, df2, col.names) {
	colnames(df2)<-col.names
	df_tmp<-rbind(df1, df2)
	return(df_tmp)
}

## --------------------
#       PARAMETERS
## --------------------
option_list = list(
		# All arguments are compulsory 
		make_option(c("-x", "--xml"), type = "character", default = NULL, metavar = "character", 
				help = "Path to the DemultiplexingStats.xml file."),
		make_option(c("-i", "--indexNumber"), type = "character", default = NULL, metavar = "character", 
				help = "Path to the .indexNumber file."),
		make_option(c("-d", "--demuxSum"), type = "character", default = NULL, metavar = "character", 
				help = "Path to the demuxSummary.txt file."),
		make_option(c("-t", "--threshold"), type = "numeric", default = 0.8, metavar = "numeric", 
				help = "Barcode count threshold")
)

opt_parser = OptionParser(usage="Make demultiplexStats easier to read.", option_list = option_list)
opt = parse_args(opt_parser)

if(is.null(opt$xml) | is.null(opt$indexNumber) | is.null(opt$demuxSum)) {
	stop("At least one argument is missing.\n", call. = FALSE)
}

## --------------------
#          LOG
## --------------------
cat("\nLancement du script demuxStatsFromXML.R avec les options suivantes :\n")
cat(paste0("\tFichier XML :\t\t", opt$xml, "\n"))
cat(paste0("\tFichier IndexNumber :\t", opt$indexNumber, "\n"))
cat(paste0("\tDemux Summary :\t\t" , opt$demuxSum, "\n"))
cat(paste0("\tBC count threshold :\t" , opt$threshold, "\n"))
launchDir<-getwd()
cat(paste0("\nLe fichier de sortie sera écrit dans le répertoire :\t",launchDir , "\n\n"))

## --------------------
#          MAIN
## --------------------
xml<-read_xml(opt$xml)

df<-data.frame()
vec.names<-c("Project", "Sample", "Barcode", "bcCount", "bcPerfect", "bcOneMismatch")

projects<-xml_find_all(xml, "//Project")

cat("Lecture du XML\n")
for (pr in 1:length(projects)){
	project<-xml_attr(projects[pr], "name")
	Samples<-xml_children(projects[pr])
	for (sample in 1:length(Samples)){
		sample_name<-xml_attr(Samples[sample], "name")
		xml_bc<-xml_children(Samples[sample])
		barcode_names<-xml_attr(xml_bc, "name")
		for (bc in 1:length(barcode_names)) {
			if (barcode_names[bc] != "all"){
				lane_path<-xml_path(xml_children(xml_bc[bc]))
				BarcodeCount<-xml_text(xml_find_all(xml, paste0(lane_path,"/BarcodeCount")))
				PerfectBarcodeCount<-xml_text(xml_find_all(xml, paste0(lane_path,"/PerfectBarcodeCount")))
				if (length(PerfectBarcodeCount) == 0) { PerfectBarcodeCount<-0 }
				OneMismatchBarcodeCount<-xml_text(xml_find_all(xml, paste0(lane_path,"/OneMismatchBarcodeCount")))
				
				if (length(OneMismatchBarcodeCount) == 0) { OneMismatchBarcodeCount<- NA}
			
				df_to_add<-data.frame(project, sample_name, barcode_names[bc], BarcodeCount, PerfectBarcodeCount, OneMismatchBarcodeCount)
				df<-concat_df(df, df_to_add, vec.names)

			}
		}		
	}
}

cat("Résumé des informations extraites (nombre d'échantillons par projet) :")
table(df$Project)

# Concaténation des index multilples
# Ecrire script pour générer ce fichier à partir de la SS
cat("\nLecture du fichier contenant le nombre d'index pour chaque échantillon.\n")
indexNumber<-read.table(opt$indexNumber, header=TRUE, sep="\t")

df2<-data.frame()
df.defaultLine<-df[which(df$Project == "default"),]
df2<-concat_df(df2, df.defaultLine, vec.names)

cat("Rassemblement des statistiques par échantillons.\n")
for (line in 1:dim(indexNumber)[1]){
	mySample<-indexNumber[line, "Sample"]
	mySampleNumber<-indexNumber[line, "NumberOfIndex"]
	cat("\nEtude de l'échantillon : " , mySample, "(" , mySampleNumber, "index )\n")
	# Single Index Case
	if (mySampleNumber == 1) {
		df.singleLine<-df[which(df$Sample == mySample),]
		df2<-concat_df(df2, df.singleLine, vec.names)
	}
	# Dual et 4 Index Cases
	else if (mySampleNumber > 1) {
		sampleName.suffixe <- ''
		if (mySampleNumber == 4) {
			sampleName.suffixe <- '_S.*'
		}
		sub.df<-df[which(str_detect(df$Sample, paste0(mySample, sampleName.suffixe))), ]
		#print(sub.df)
		if (nrow(sub.df) == 0) {
			cat("Aucun échantillon trouvé !\n")
			cat("La recherche de l'échantillon", paste0(mySample, sampleName.suffixe), "dans le data.table suivant à échouée :\n")
			print(df)
		} else {
			countBarcodesDone = 1
			# Parcours du sous-data.frame
			for (l in 1:dim(sub.df)[1]) {
				sub.df.project<-sub.df[l, "Project"]
				sub.df.barcode<-sub.df[l, "Barcode"]
				sub.df.bcCount<-as.numeric(sub.df[l, "bcCount"])
				sub.df.bcPerfect<-as.numeric(sub.df[l, "bcPerfect"])
				sub.df.oneMismatch<-as.numeric(sub.df[l, "bcOneMismatch"])	# bcOneMismatch
				
				# Première iteration
				countBarcodesDone = countBarcodesDone + str_count(sub.df.barcode, "\\+")
				if (countBarcodesDone <= mySampleNumber) {
					if (l == 1 ) {
						sub.df.project.toAdd<-sub.df.project
						sub.df.barcode.toAdd<-sub.df.barcode
						sub.df.bcCount.toAdd<-sub.df.bcCount
						sub.df.bcPerfect.toAdd<-sub.df.bcPerfect
						sub.df.oneMismatch.toAdd<-sub.df.oneMismatch
					} else {
						sub.df.barcode.toAdd<-paste0(sub.df.barcode.toAdd, "+", sub.df.barcode)
						sub.df.bcCount.toAdd<-sub.df.bcCount.toAdd+sub.df.bcCount
						sub.df.bcPerfect.toAdd<-sub.df.bcPerfect.toAdd+sub.df.bcPerfect
						sub.df.oneMismatch.toAdd<-sub.df.oneMismatch.toAdd+sub.df.oneMismatch
					}
				}
			}
			# Add to data.frame
			df_to_add<-data.frame(sub.df.project,mySample, sub.df.barcode.toAdd, sub.df.bcCount.toAdd, sub.df.bcPerfect.toAdd, sub.df.oneMismatch.toAdd)
			df2<-concat_df(df2, df_to_add, vec.names)
		}	
	}
}
	
cat("\nRésumé des informations extraites (nombre d'échantillons par projet) :")
table(df2$Project)

## Recherche des index indeterminés
cat("\nRecherche des index indéterminés.\n")
bcCount.min<-min(as.numeric(df2[-which(df$Project == "default"), "bcCount"]))
bcCount.threshold<-opt$threshold*bcCount.min

# Rechercher tous les index trouvés au moins bcCount.threshold fois
cat("Tentative de récupérer des échantillons parmi les index retrouvés les plus fréquemment.\n")
cat("\tLecture du DemuxSummary.\n")
linesToSkip<-as.numeric(system(paste("grep -n Most", opt$demuxSum, "| cut -d':' -f1"), intern = TRUE))
tabDemuxSum<-read.table(opt$demuxSum, skip=linesToSkip, col.names=c("Index", "Count"))

tabUndetermined<-tabDemuxSum[which(tabDemuxSum$Count >= bcCount.threshold),]

cat("\tRésumé des inforamtions extraites :\n")
cat(paste0("\tNombre d'index indéterminés retrouvés :\t", dim(tabUndetermined)[1], "\n"))
if(nrow(tabUndetermined) > 0) { head(tabUndetermined) }


# Construction du dataFrame pour intégration à df2
#df2.Projects<-unique(df2$Project)
#myProject<-df2.Projects[which(df2.Projects != "default")]

### Pour chaque ligne de tabUndertermined, on ajoute une ligne à df2 :
if (dim(tabUndetermined)[1] != 0) {
	df.tabUndetermined<-data.frame()
	for (i in 1:dim(tabUndetermined)[1]) {
		df.tabUndetermined.tmp<-data.frame("default", "Undetermined", tabUndetermined[i, "Index"], tabUndetermined[i, "Count"], NA, NA)
		df.tabUndetermined<-concat_df(df.tabUndetermined, df.tabUndetermined.tmp, vec.names)
	}
	
	df2<-concat_df(df2, df.tabUndetermined, vec.names)
	cat("\tLes index indéterminés ont été ajouté au data.table.\n")
} else {
	cat("\tAuncun index indéterminés trouvés.\n")
}

## Soustraction des undertermined aux allOthers
# recuperer les Count de tabUndetermined et soustraire la somme à df2[which(df2$Barcode == "unknown"), "bcCount"]
cat("\nQuelques calculs sur les données avant de les exporter.\n")
cat("\tActualisation du nombre d'index 'AllOthers'.\n")
undertermined.count<-sum(as.numeric(tabUndetermined[,"Count"]))
df2[which(df2$Barcode == "unknown"), "bcCount"]<-as.numeric(df2[which(df2$Barcode == "unknown"), "bcCount"])-undertermined.count

# Calcul pourcentages de chaque barcode
cat("\tCalcul du pourcentage sur le nombre de fragments total.\n")
totalOfFragments<-sum(as.numeric(df2$bcCount))

percentOfFragment<-as.data.frame(round((as.numeric(df2[,"bcCount"])/totalOfFragments)*100, 2))
rownames(percentOfFragment)<-rownames(df2)
colnames(percentOfFragment)<-"percentageOfFragment"

df2<-cbind(df2, percentOfFragment)

# Remplacement des NA par 0
df2[is.na(df2)] <- 0

# Export du data.frame
cat("\nSauvegarde du data.frame.\n")
#myProject<-"DEBUG"
write.table(df2, row.names = FALSE, quote = F, sep = "\t", file = paste0("DemultiplexStats.tsv"))
# Ecrire un fichier par valeur de myProject ! Cas ou il y a plusieurs projets sur la même lane.
cat(paste0("\tLe fichier suivant à été créé :\t", launchDir, "/DemultiplexStats.tsv\n"))
cat("\nFin normale du script, on sort.\n")