Wrote by Huamei Li
State Key Laboratory of Bioelectronics, School of Biological Science & Medical Engineering, Southeast University, Nanjing 210096, China;
Source: DeconPeaker_demo.ipynb
Here, we will demonstrate how DeconPeaker works on chromatin accessibility datasets (ATAC-Seq) as well as gene expression datasets (RNA-Seq & Microarray), and employing Cores'syntentic dataset (GSE74912, ATAC-Seq), Liu's dataset (GSE64098, RNS-Seq) and Shen-Orr's dataset (GSE49830, Microarray) as examples.
DeconPeaker provides a commond line to process chromatin accessibility datasets except gene expression datasets, and this step only supports Linux system. The input file must be YAML format, and including five keys. 1) sample_ID: unique sample ID for each sample.2) bam: BAM of aligned file path 3) batch: batch ID of each sample, if cell type has replicates, batch ID must be the same.4) cellname: cell name of the pure sample. *5) peak: peak file path of the sample. The test YAML format file can see in test/GSE74912_ATAC_pure_cells.yaml
# python deconPeaker.py preprocess --pure=test/GSE74912_ATAC_pure_cells.yaml --thread=10 --prefix=Test --outdir=results
This commond line will create a directory, named "preprocess", including count matrix file and phenotype class file. Both files are the basis for the subsequent analysis.
import matplotlib
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns; sns.set()
import matplotlib.image as mpimg
%matplotlib inline
# Signature matrix file of ATAC-Seq data
profile = pd.read_csv("./test/examples/ATAC-Seq/GSE74912_Corces_MR_pure_readcounts_signature_matrix.xls", sep='\t', header=0)
# first column: chromosome name; second column: start position of the peak; third column: end position of the peak
profile.head()
mixtures = pd.read_csv('test\examples\ATAC-Seq\GSE74912_Corces_MR_synthetic_mixture_counts_data.xls', sep='\t', header=0)
mixtures.head()
# Estimate P value of the deconvolution
%run -i deconPeaker.py deconvolution --lib-strategy=ATAC-Seq --mixture=test\examples\ATAC-Seq\GSE74912_Corces_MR_synthetic_mixture_counts_data.xls --pure=test\examples\ATAC-Seq\GSE74912_Corces_MR_pure_readcounts_signature_matrix.xls --format=TABLE --pvalue=FALSE --outdir=results\GSE74912_Corces_MR
%run -i test\true_VS_est_scatter_plot.py test\examples\ATAC-Seq\GSE74912_Corces_MR_Dirichlet_proportions.xls results\GSE74912_Corces_MR\deconvolution\deconPeaker-Results.xls results\GSE74912_Corces_MR\deconvolution\Truths_Vs_Predictions
# input file of findstsps step of DeconPeaker, of which rows represent genes and columns sample names.
profile = pd.read_csv("./test/examples/RNA-Seq/GSE64098_liu_phenotype_classes.xls", sep='\t', header=0, index_col=0)
profile.head()
phenotype = pd.read_csv("./test/examples/RNA-Seq/GSE64098_liu_phenotype_classes.xls", sep='\t', header=0, index_col=0)
# value 1 = indicate membership of the reference sample, value 2 = indicates the class that the sample will be compared against.
phenotype.head()
%run -i deconPeaker.py findctsps --lib-strategy=RNA-Seq --profile=./test/examples/RNA-Seq/GSE64098_liu_pure_counts_data.xls --phenotype=./test/examples/RNA-Seq/GSE64098_liu_phenotype_classes.xls --outdir=results/GSE64098_liu
signature_matrix = pd.read_csv('results/GSE64098_liu/findctsps/GSE64098_liu_pure_counts_data_signature_matrix.xls', sep='\t', index_col=0)
signature_matrix.head()
img = mpimg.imread('results/GSE64098_liu/findctsps/GSE64098_liu_pure_counts_data_signature_heatmap.png')
plt.imshow(img)
plt.show()
%run -i deconPeaker.py deconvolution --lib-strategy=RNA-Seq --mixture=test\examples\RNA-Seq\GSE64098_liu_mixture_counts_data.xls --pure=results\GSE64098_liu\findctsps\GSE64098_liu_pure_counts_data_signature_matrix.xls --format=TABLE --pvalue=TRUE --outdir=results\GSE64098_liu
%run -i test\true_VS_est_scatter_plot.py test\examples\RNA-Seq\GSE64098_liu_annotation_mixture.xls results\GSE64098_liu\deconvolution\deconPeaker-Results.xls results\GSE64098_liu\deconvolution\Truths_Vs_Predictions
%run -i deconPeaker.py findctsps --lib-strategy=Microarray --profile=./test/examples/Microarray/GSE19830_Shen_Orr_pure_data.xls --phenotype=./test/examples/Microarray/GSE19830_Shen_Orr_phenotype_classes.xls --outdir=results/GSE19830_Shen_Orr
signature_matrix = pd.read_csv('results/GSE19830_Shen_Orr/findctsps/GSE19830_Shen_Orr_pure_data_signature_matrix.xls', sep='\t', index_col=0)
signature_matrix.head()
img = mpimg.imread('results/GSE19830_Shen_Orr/findctsps/GSE19830_Shen_Orr_pure_data_signature_heatmap.png')
plt.imshow(img)
plt.show()
# Run DeconPeaker without estimating P-value
%run -i deconPeaker.py deconvolution --lib-strategy=Microarray --mixture=test\examples\Microarray\GSE19830_Shen_Orr_mixture_data.xls --pure=results\GSE19830_Shen_Orr\findctsps\GSE19830_Shen_Orr_pure_data_signature_matrix.xls --format=TABLE --pvalue=FALSE --outdir=results\GSE19830_Shen_Orr
%run -i test\true_VS_est_scatter_plot.py test\examples\Microarray\GSE19830_Shen_Orr_annotation_mixture.xls results\GSE19830_Shen_Orr\deconvolution\deconPeaker-Results.xls results\GSE19830_Shen_Orr\deconvolution\Truths_Vs_Predictions