Merge pull request #24 from luslab/dev (48911644) · Commits · Chaos / Cutandrun

bin/consensus_peaks.py

0 → 100755

+74 −0

Original line number	Diff line number	Diff line
		#!/usr/bin/env python

		import os
		import glob
		import re
		import numpy as np
		import pandas as pd
		import matplotlib.pyplot as plt
		import argparse
		import upsetplot

		############################################
		############################################
		## PARSE ARGUMENTS
		############################################
		############################################
		Description = 'Upset ven diagram of consensus peaks.'
		Epilog = """Example usage: python consensus_peaks.py <MERGED_INTERVAL_FILE> """

		parser = argparse.ArgumentParser(description=Description, epilog=Epilog)

		## REQUIRED PARAMETERS
		parser.add_argument('--peaks', help="Merged peaks interval file with replicate counts column.")
		parser.add_argument('--outpath', help="Full path to output directory.")
		args = parser.parse_args()

		############################################
		############################################
		## MAIN FUNCTION
		############################################
		############################################

		# create list of data frames, one for each group consensus peaks file
		peak_file_list = glob.glob(args.peaks)
		peak_df_list = list()
		for i in list(range(len(peak_file_list))):
		peaks_i = pd.read_csv(peak_file_list[i], sep='\t', header=None, usecols=[0,1,2,8,9], names=['chrom','start','end','sample_reps','count'])
		peaks_i['sample_reps'] = peaks_i['sample_reps'].replace(".peaks.bed.stringent.bed", "", regex=True)
		peak_df_list.append(peaks_i)
		reps2 = peaks_i[peaks_i["count"]>1]

		# add sorted column to each dataframe, and make new condensed dataframe list
		summary_peak_df_list = list()
		for i in list(range(len(peak_df_list))):
		peaks_i = peak_df_list[i]
		peaks_i['sorted_samples'] = ''
		rows_now = peaks_i.shape[0]
		for j in list(range(rows_now)):
		sample_list = peaks_i.at[j,'sample_reps']
		sample_array = np.unique(sample_list.split(','))
		sample_sorted = sorted(sample_array)
		sample_str = ",".join(sample_sorted)
		peaks_i.at[j,'sorted_samples'] = sample_str
		summary_peaks_i = peaks_i[['sorted_samples', 'count']].groupby(['sorted_samples'], as_index = False).sum()
		summary_peak_df_list.append(summary_peaks_i)

		# construct data in appropriate format for upsetplot, and plot
		for i in list(range(len(summary_peak_df_list))):
		df_i = summary_peak_df_list[i]
		# Get group name
		basename = os.path.basename(peak_file_list[i])
		group_name = basename.rsplit(".", -1)[0]
		file_name = group_name + ".consensus_peaks.pdf"
		categories = df_i.shape[0]
		cat_list = []
		for j in list(range(categories)):
		summary_sample = df_i.at[j,'sorted_samples'].split(',')
		cat_list.append(summary_sample)

		# Plot
		peak_counts = upsetplot.from_memberships(cat_list, data = df_i['count'])
		upsetplot.plot(peak_counts)
		plt.show()
		plt.savefig(os.path.join(args.outpath, file_name))

bin/deseq2_diff.r

+41 −184

Original line number	Diff line number	Diff line
		@@ -57,16 +57,6 @@ if (is.null(opt$groups)){
		stop("Please provide group names for analysis.", call.=FALSE)
		}

		#if (is.null(opt$control)){
		# print_help(opt_parser)
		# stop("Please provide a control group name.", call.=FALSE)
		#}
		#
		#if (is.null(opt$treatment)){
		# print_help(opt_parser)
		# stop("Please provide a treatment group name.", call.=FALSE)
		#}

		if (is.null(opt$bed)){
		print_help(opt_parser)
		stop("Please provide a list of bed files to load.", call.=FALSE)
		@@ -81,12 +71,6 @@ if (is.null(opt$bam)){
		bed_list <- unlist(strsplit(opt$bed, ","))
		bam_list <- unlist(strsplit(opt$bam, ","))

		# Check same length <- MAYBE TAKE THIS OUT IF WE SIMPLY PASS THE ALIGNMENT CHANNEL AS IT WILL CONTAIN IGG
		#if (length(bed_list) != length(bam_list)) {
		# print_help(opt_parser)
		# stop("Bed and bam file list are different lengths", call.=FALSE)
		#}

		################################################
		################################################
		## READ IN DATA FILES ##
		@@ -105,10 +89,18 @@ if (!is.null(opt$exclude)) {
		groups = groups[!matching]
		}

		if (length(unique(groups)) == 1) {
		message("WARN: Only one experimental group identified, so differential analysis cannot be executed. DESEQ2 will be skipped.")
		quit(save = "no", status = 0, runLast = FALSE)
		}

		################ redo from here #######################

		# Init
		mPeak = GRanges()
		file_count = 0
		file_list=vector()
		sample_mat = matrix(NA, length(bam_list), 2)
		# Read in bed files that match the control or treatment group
		for(group in groups){
		search_res <- str_detect(bed_list, group)
		@@ -117,39 +109,40 @@ for(group in groups){
		}
		file_list <- bed_list[search_res] %>% append(file_list, .)
		}
		k=1
		for(file in file_list) {
		peakRes = read.table(file, header = FALSE, fill = TRUE)
		mPeak = GRanges(seqnames = peakRes$V1, IRanges(start = peakRes$V2, end = peakRes$V3), strand = "*") %>% append(mPeak, .)
		file_count = file_count + 1
		file_name = basename(file)
		file_name = strsplit(file_name, "[.]")[[1]][1]
		file_split = strsplit(file_name, "[_]")[[1]]
		file_group = paste(file_split[-length(file_split)], collapse = "_")
		file_rep = file_split[length(file_split)]
		file_rep_num = gsub("[^0-9.]", "", file_rep)
		sample_mat[k,1] = file_group
		sample_mat[k,2] = file_rep_num
		k=k+1
		}



		# Create replicate counts and names
		group_count = length(groups)
		rep_count = file_count / group_count
		reps = paste0("rep", 1:rep_count)
		expected_rep_str = paste0("_R", 1:rep_count)

		# Create peak table and count matrix
		masterPeak = reduce(mPeak)
		countMat = matrix(NA, length(masterPeak), file_count)
		colnames(countMat) = paste(rep(groups, rep_count), rep(reps, each = group_count, sep = "_"))

		sample_list = paste(sample_mat[,1], sample_mat[,2], sep = "_R")
		colnames(countMat) = sample_list
		# Read in bam files that match the control or treatment group
		for(i in seq_along(groups)){
		search_res <- str_detect(bam_list, groups[i])
		if (!any(search_res)) {
		stop(paste("group", i, "was not found amongst bam files"))
		}
		file_list <- bam_list[search_res]
		}

		for(j in seq_along(file_list)) {
		rep_search <- str_detect(file_list, expected_rep_str[j])
		file_now <- file_list[rep_search]
		for(j in seq_along(bam_list)) {
		rep_search <- str_detect(bam_list, sample_list[j])
		file_now <- bam_list[rep_search]
		fragment_counts <- getCounts(file_now, masterPeak, paired = TRUE, by_rg = FALSE, format = "bam")
		countMat[, ((j*group_count) - (group_count - i))] = counts(fragment_counts)[,1]
		}
		countMat[,j] = counts(fragment_counts)[,1]
		}

		################################################
		@@ -162,20 +155,28 @@ if (file.exists(opt$outdir) == FALSE) {
		dir.create(opt$outdir,recursive=TRUE)
		}
		setwd(opt$outdir)
		#if (FALSE) {

		## Create index list for peak count filter
		selectR = which(rowSums(countMat) > opt$count_thresh)
		dataS = countMat[selectR,] ## Select data from filter
		condition = factor(rep(groups, each = length(reps)))
		condition = factor(sample_mat[,1])

		samples.vec <- sort(colnames(countMat))
		#groups <- sub("_[^_]+$", "", samples.vec)
		if (length(unique(groups)) == 1 \|\| length(unique(groups)) == length(samples.vec)) {

		if (length(unique(groups)) == length(samples.vec)) {
		message("WARN: Only one experimental replicate per group was identified, so differential analysis cannot be executed. DESEQ2 will be skipped.")
		quit(save = "no", status = 0, runLast = FALSE)
		}

		DDSFile <- paste(opt$outprefix,".dds.RData",sep="")
		if (file.exists(DDSFile) == FALSE) {
		tryCatch (
		expr = DESeqDataSetFromMatrix(countData = dataS,colData = DataFrame(condition),design = ~ condition),
		error = function(cond){
		message("failed to create DESeqDataSet object. DESEQ2 module will be skipped")
		quit(save = "no", status = 0, runLast = FALSE)
		}
		)
		dds = DESeqDataSetFromMatrix(countData = dataS,colData = DataFrame(condition),design = ~ condition)
		dds <- DESeq(dds, parallel=TRUE, BPPARAM=MulticoreParam(opt$cores))
		if (!opt$vst) {
		@@ -206,49 +207,6 @@ res = results(DDS, independentFiltering = FALSE, altHypothesis = "greaterAbs")
		countMatDiff = cbind(dataS, normDDS, res) ## Combine deseq results
		write.csv(countMatDiff, paste(opt$outprefix, "count_mat_diff.csv", sep=""), row.names=FALSE)


		################################################
		################################################
		## NEW PLOT QC ##
		################################################
		################################################
		if (FALSE) {
		PlotFile <- paste(opt$outprefix,".plots.pdf",sep="")
		if (file.exists(PlotFile) == FALSE) {
		pdf(file=PlotFile,onefile=TRUE,width=7,height=7)

		## PCA
		pca.data <- DESeq2::plotPCA(rld,intgroup=c("condition"),returnData=TRUE)
		percentVar <- round(100 * attr(pca.data, "percentVar"))
		plot <- ggplot(pca.data, aes(PC1, PC2, color=condition)) +
		geom_point(size=3) +
		xlab(paste0("PC1: ",percentVar[1],"% variance")) +
		ylab(paste0("PC2: ",percentVar[2],"% variance")) +
		theme(panel.grid.major = element_blank(),
		panel.grid.minor = element_blank(),
		panel.background = element_blank(),
		panel.border = element_rect(colour = "black", fill=NA, size=1))
		print(plot)

		## WRITE PC1 vs PC2 VALUES TO FILE
		pca.vals <- pca.data[,1:2]
		colnames(pca.vals) <- paste(colnames(pca.vals),paste(percentVar,'% variance',sep=""), sep=": ")
		pca.vals <- cbind(sample = rownames(pca.vals), pca.vals)
		write.table(pca.vals,file=paste(opt$outprefix,".pca.vals.txt",sep=""),row.names=FALSE,col.names=TRUE,sep="\t",quote=TRUE)

		## SAMPLE CORRELATION HEATMAP
		sampleDists <- dist(t(assay(rld)))
		sampleDistMatrix <- as.matrix(sampleDists)
		colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
		pheatmap(sampleDistMatrix,clustering_distance_rows=sampleDists,clustering_distance_cols=sampleDists,col=colors)

		## WRITE SAMPLE DISTANCES TO FILE
		write.table(cbind(sample = rownames(sampleDistMatrix), sampleDistMatrix),file=paste(opt$outprefix,".sample.dists.txt",sep=""),row.names=FALSE,col.names=TRUE,sep="\t",quote=FALSE)

		dev.off()
		}
		}

		################################################
		################################################
		## PLOT QC ##
		@@ -305,6 +263,7 @@ if (file.exists(PlotFile) == FALSE) {
		pca.data <- plotPCA_vst(dds, assay=vst_name,intgroup="condition",ntop=n_top_var)
		percentVar <- round(attr(pca.data, "percentVar")$percentVar)
		plot_subtitle <- ifelse(n_top_var==Inf, "All peaks", paste("Top", n_top_var, "peaks"))
		# PCA PLOT 1 - FIRST PC
		pl <- ggplot(pca.data, aes(PC1, PC2, color=condition)) +
		geom_point(size=3) +
		xlab(paste0("PC1: ",percentVar[1],"% variance")) +
		@@ -317,6 +276,7 @@ if (file.exists(PlotFile) == FALSE) {
		panel.border = element_rect(colour = "black", fill=NA, size=1))
		print(pl)

		# PCA PLOT 2 - DIAGNOSTIC OF PCs
		pl <- ggplot(attr(pca.data, "percentVar"), aes(x=PC, y=percentVar)) +
		geom_line(aes(colour="explained by PC")) +
		geom_line(aes(y=groupR, colour="of PC explained by condition")) +
		@@ -326,6 +286,7 @@ if (file.exists(PlotFile) == FALSE) {
		theme(legend.position="top")
		print(pl)

		# PCA PLOT 3 - GROUP-EXPLANATORY PCs
		pc_r <- order(attr(pca.data, "percentVar")$groupR, decreasing=TRUE)
		pl <- ggplot(pca.data, aes_string(paste0("PC", pc_r[1]), paste0("PC", pc_r[2]), color="condition")) +
		geom_point(size=3) +
		@@ -367,109 +328,6 @@ if (file.exists(PlotFile) == FALSE) {
		dev.off()
		}

		################################################
		################################################
		## LOOP THROUGH COMPARISONS ##
		################################################
		################################################

		if (FALSE) {
		ResultsFile <- paste(opt$outprefix,".results.txt",sep="")
		if (file.exists(ResultsFile) == FALSE) {
		raw.counts <- counts(dds,normalized=FALSE)
		colnames(raw.counts) <- paste(colnames(raw.counts),'raw',sep='.')
		pseudo.counts <- counts(dds,normalized=TRUE)
		colnames(pseudo.counts) <- paste(colnames(pseudo.counts),'pseudo',sep='.')

		deseq2_results_list <- list()
		comparisons <- combn(unique(groups),2)
		for (idx in 1:ncol(comparisons)) {
		control.group <- comparisons[1,idx]
		treat.group <- comparisons[2,idx]
		CompPrefix <- paste(control.group,treat.group,sep="vs")
		cat("Saving results for ",CompPrefix," ...\n",sep="")

		CompOutDir <- paste(CompPrefix,'/',sep="")
		if (file.exists(CompOutDir) == FALSE) {
		dir.create(CompOutDir,recursive=TRUE)
		}

		control.samples <- samples.vec[which(groups == control.group)]
		treat.samples <- samples.vec[which(groups == treat.group)]
		comp.samples <- c(control.samples,treat.samples)

		comp.results <- results(dds,contrast=c("condition",c(control.group,treat.group)))
		comp.df <- as.data.frame(comp.results)
		#comp.table <- cbind(interval.table, as.data.frame(comp.df), raw.counts[,paste(comp.samples,'raw',sep='.')], pseudo.counts[,paste(comp.samples,'pseudo',sep='.')])
		comp.table <- cbind(as.data.frame(comp.df), raw.counts[,paste(comp.samples,'raw',sep='.')], pseudo.counts[,paste(comp.samples,'pseudo',sep='.')])

		## WRITE RESULTS FILE
		CompResultsFile <- paste(CompOutDir,CompPrefix,opt$outsuffix,".deseq2.results.txt",sep="")
		write.table(comp.table, file=CompResultsFile, col.names=TRUE, row.names=FALSE, sep='\t', quote=FALSE)

		## FILTER RESULTS BY FDR & LOGFC AND WRITE RESULTS FILE
		pdf(file=paste(CompOutDir,CompPrefix,opt$outsuffix,".deseq2.plots.pdf",sep=""),width=10,height=8)
		if (length(comp.samples) > 2) {
		for (MIN_FDR in c(0.01,0.05)) {

		## SUBSET RESULTS BY FDR
		#pass.fdr.table <- subset(comp.table, padj < MIN_FDR)
		#pass.fdr.up.table <- subset(comp.table, padj < MIN_FDR & log2FoldChange > 0)
		#pass.fdr.down.table <- subset(comp.table, padj < MIN_FDR & log2FoldChange < 0)

		## SUBSET RESULTS BY FDR AND LOGFC
		#pass.fdr.logFC.table <- subset(comp.table, padj < MIN_FDR & abs(log2FoldChange) >= 1)
		#pass.fdr.logFC.up.table <- subset(comp.table, padj < MIN_FDR & abs(log2FoldChange) >= 1 & log2FoldChange > 0)
		#pass.fdr.logFC.down.table <- subset(comp.table, padj < MIN_FDR & abs(log2FoldChange) >= 1 & log2FoldChange < 0)

		## WRITE RESULTS FILE
		#CompResultsFile <- paste(CompOutDir,CompPrefix,opt$outsuffix,".deseq2.FDR",MIN_FDR,".results.txt",sep="")
		#CompBEDFile <- paste(CompOutDir,CompPrefix,opt$outsuffix,".deseq2.FDR",MIN_FDR,".results.bed",sep="")
		#write.table(pass.fdr.table, file=CompResultsFile, col.names=TRUE, row.names=FALSE, sep='\t', quote=FALSE)
		#write.table(pass.fdr.table[,c("Chr","Start","End","Geneid","log2FoldChange","Strand")], file=CompBEDFile, col.names=FALSE, row.names=FALSE, sep='\t', quote=FALSE)

		## MA PLOT & VOLCANO PLOT
		DESeq2::plotMA(comp.results, main=paste("MA plot FDR <= ",MIN_FDR,sep=""), ylim=c(-2,2),alpha=MIN_FDR)
		#plot(comp.table$log2FoldChange, -1log10(comp.table$padj), col=ifelse(comp.table$padj<=MIN_FDR, "red", "black"), xlab="logFC", ylab="-1log10(FDR)", main=paste("Volcano plot FDR <=",MIN_FDR,sep=" "), pch=20)

		## ADD COUNTS TO LOGFILE
		#cat(CompPrefix," genes with FDR <= ",MIN_FDR,": ",nrow(pass.fdr.table)," (up=",nrow(pass.fdr.up.table),", down=",nrow(pass.fdr.down.table),")","\n",file=LogFile,append=TRUE,sep="")
		#cat(CompPrefix," genes with FDR <= ",MIN_FDR," & FC > 2: ",nrow(pass.fdr.logFC.table)," (up=",nrow(pass.fdr.logFC.up.table),", down=",nrow(pass.fdr.logFC.down.table),")","\n",file=LogFile,append=TRUE,sep="")

		}
		cat("\n",file=LogFile,append=TRUE,sep="")
		}

		## SAMPLE CORRELATION HEATMAP
		rld.subset <- assay(rld)[,comp.samples]
		sampleDists <- dist(t(rld.subset))
		sampleDistMatrix <- as.matrix(sampleDists)
		colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
		pheatmap(sampleDistMatrix,clustering_distance_rows=sampleDists,clustering_distance_cols=sampleDists,col=colors)

		## SCATTER PLOT FOR RLOG COUNTS
		combs <- combn(comp.samples,2,simplify=FALSE)
		clabels <- sapply(combs,function(x){paste(x,collapse=' & ')})
		plotdat <- data.frame(x=unlist(lapply(combs, function(x){rld.subset[, x[1] ]})),y=unlist(lapply(combs, function(y){rld.subset[, y[2] ]})),comp=rep(clabels, each=nrow(rld.subset)))
		plot <- xyplot(y~x\|comp,plotdat,
		panel=function(...){
		panel.xyplot(...)
		panel.abline(0,1,col="red")
		},
		par.strip.text=list(cex=0.5))
		print(plot)
		dev.off()

		colnames(comp.df) <- paste(CompPrefix,".",colnames(comp.df),sep="")
		deseq2_results_list[[idx]] <- comp.df
		}

		## WRITE RESULTS FROM ALL COMPARISONS TO FILE
		deseq2_results_table <- cbind(do.call(cbind, deseq2_results_list),raw.counts,pseudo.counts)
		write.table(deseq2_results_table, file=ResultsFile, col.names=TRUE, row.names=FALSE, sep='\t', quote=FALSE)
		}
		}

		################################################
		################################################
		## SAVE SIZE FACTORS ##
		@@ -507,7 +365,6 @@ if (file.exists(RLogFile) == FALSE) {
		print(a)
		sink()
		}
		#}
		################################################
		################################################
		################################################

conf/modules.config

+24 −3

Original line number	Diff line number	Diff line
		@@ -41,6 +41,12 @@ params {
		"deseq2" {
		publish_dir = "deseq2_qc"
		}
		"bedtools_merge_groups" {
		args = " -c 2,3,4,5,6,7,7 -o collapse,collapse,collapse,collapse,collapse,collapse,count_distinct"
		// args = " -c 7 -o collapse"
		publish_dir = "seacr/consensus_peaks"
		suffix = ".consensus.peaks"
		}
		"bedtools_genomecov_bedgraph" {
		args = ""
		suffix = ""
		@@ -51,6 +57,10 @@ params {
		suffix = ".peaks.bed"
		publish_dir = "seacr"
		}
		"sort_group_peaks" {
		args = "-k1,1 -k2,2n"
		publish_files = false
		}
		"ucsc_bedclip" {
		suffix = ".clip"
		publish_files = false
		@@ -134,6 +144,9 @@ params {
		args = "--sortUsing sum --startLabel \"Peak Start\" --endLabel \"Peak End\" --xAxisLabel \"\" --regionsLabel \"Peaks\""
		publish_dir = "deeptools/heatmaps/peaks"
		}
		"plot_peaks" {
		publish_dir = "seacr"
		}

		/*
		========================================================================================
		@@ -181,6 +194,14 @@ params {
		suffix = ".max_signal"
		publish_files = false
		}
		"awk_name_peak_bed" {
		command = "'{OFS = \"\\t\"} {print \$0, FILENAME}'"
		publish_files = false
		}
		"awk_threshold" {
		publish_dir = "seacr/consensus_peaks"
		suffix = ".rep_thresh"
		}

		/*
		========================================================================================

dev/docker/consensus_peaks/Dockerfile

0 → 100644

+8 −0

Original line number	Diff line number	Diff line
		FROM nfcore/base:1.14
		LABEL authors="charlotte.west@crick.ac.uk" \
		description="Docker image containing all requirements for development of thresholding peaks and plotting"

		# Install conda packages
		COPY ./environment.yml /
		RUN conda env create -f /environment.yml && conda clean -a
		ENV PATH /opt/conda/envs/reporting/bin:$PATH

dev/docker/consensus_peaks/environment.yml

0 → 100644

+14 −0

Original line number	Diff line number	Diff line
		# conda env create -f environment.yml
		name: reporting
		channels:
		- conda-forge
		- bioconda
		- defaults
		dependencies:
		# python version
		- python=3.8.3

		# conda packages
		- numpy=1.20.*
		- pandas=1.2.*
		- upsetplot=0.4.4

Admin message