Skip to content
Snippets Groups Projects
Commit 142e5028 authored by Gervaise Henry's avatar Gervaise Henry :cowboy:
Browse files

Merge remote-tracking branch 'origin/develop' into PrFx

parents 9444c3d9 ca6ba2ef
No related merge requests found
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1478 additions and 261 deletions
.gitignore
+ 1
− 0
View file @ 142e5028
......@@ -10,3 +10,4 @@ WR/
*~
temp_png.png
Rplots.pdf
stress/*
README.md
+ 7
− 3
View file @ 142e5028
......@@ -11,8 +11,8 @@ Determining cellular heterogeneity in the human prostate with single-cell RNA se
* <a href="https://orcid.org/0000-0002-0746-927X" target="orcid.widget" rel="noopener noreferrer" style="vertical-align:top;"><img src="https://orcid.org/sites/default/files/images/orcid_16x16.png" style="width:1em;margin-right:.5em;" alt="ORCID iD icon">orcid.org/0000-0002-0746-927X</a>
* PI Email: [douglas.strand@utsouthwestern.edu](mailto:douglas.strand@utsouthwestern.edu)
* **ANALYZED DATA FOR QUERYING AT: [StrandLab.net](http://strandlab.net/analysis.php)**
* **Raw data at: [GEO (scRNA-Seq)](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117403) & [GEO (popRNA-Seq)](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117271) & [GenitoUrinary Development Molecular Anatomy Project (GUDMAP)]("https://doi.org/10.25548/W-R8CM")**
* **Publication at: PENDING**
* **Raw data at: [GEO](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE120716) & [GenitoUrinary Development Molecular Anatomy Project (GUDMAP)]("https://doi.org/10.25548/W-R8CM")**
* **Publication at: [BioRxiv](https://www.biorxiv.org/content/early/2018/10/15/439935)**
Data Analysis
-------------
......@@ -50,6 +50,10 @@ Data Analysis
* run bash script [sc\_TissueMapper\-D17\_FACS.sh](https://git.biohpc.swmed.edu/StrandLab/sc-TissueMapper_Pr/blob/master/bash.scripts/sc_TissueMapper-D17_FACS.sh)
* 3 Run on 2nd patient FACS samples
* run bash script [sc\_TissueMapper\-D27\_FACS.sh](https://git.biohpc.swmed.edu/StrandLab/sc-TissueMapper_Pr/blob/master/bash.scripts/sc_TissueMapper-D27_FACS.sh)
* 4 Run on several downsamples from 1 sample from 1st patient
* run bash script [sc\_TissueMapper\-DS\_D17.sh](https://git.biohpc.swmed.edu/StrandLab/sc-TissueMapper_Pr/blob/master/bash.scripts/sc_TissueMapper-DS_D17.sh)
* 5 Aggregate and compare several downsamples from #4
* run bash script [sc\_TissueMapper\_RUN.DS\_D17.aggr.R](https://git.biohpc.swmed.edu/StrandLab/sc-TissueMapper_Pr/blob/master/r.scripts/sc_TissueMapper_RUN.DS_D17.aggr.R)
* **Pipeline:**
* Link cellranger count/aggr output to analysis
* Create demultiplex file to add custom sample groups
......@@ -118,4 +122,4 @@ Data Analysis
* "c2.all.v6.1.symbols.gmt" MSigDB C2 Curated Gene Sets [**MSigDB C2**](http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection=C2)
* "c2.cp.kegg.v6.1.symbols" MSigDB C2 KEGG Gene Subsets [**KEGG**](http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection=CP:KEGG)
* "c5.all.v6.1.symbols.gmt" MSigDB C5 Gene Ontology Gene Sets [**MSigDB C5**](http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection=C5)
* "c5.bp.v6.1.symbols.gmt" MSigDB C5 Gene Ontology Biological Processes Gene Subsets [**GO BP**](http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection=BP)
\ No newline at end of file
* "c5.bp.v6.1.symbols.gmt" MSigDB C5 Gene Ontology Biological Processes Gene Subsets [**GO BP**](http://software.broadinstitute.org/gsea/msigdb/genesets.jsp?collection=BP)
bash.scripts/sc_TissueMapper-D17_FACS.sh
+ 4
− 2
View file @ 142e5028
#!/bin/bash
#SBATCH --job-name R_FullAnalysis
#SBATCH --job-name R_FullAnalysis.D17FACS
#SBATCH -p 256GB,256GBv1,384GB
#SBATCH -N 1
#SBATCH -t 7-0:0:0
......@@ -8,6 +8,8 @@
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load R/3.4.1-gccmkl
module load R/3.4.1-gccmkl_20181025
sh ./sc_LinkData.sh D17_FACS
Rscript ../r.scripts/sc-TissueMapper_RUN.D17_FACS.R
bash.scripts/sc_TissueMapper-D27_FACS.sh
+ 4
− 2
View file @ 142e5028
#!/bin/bash
#SBATCH --job-name R_FullAnalysis
#SBATCH --job-name R_FullAnalysis.D27FACS
#SBATCH -p 256GB,256GBv1,384GB
#SBATCH -N 1
#SBATCH -t 7-0:0:0
......@@ -8,6 +8,8 @@
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load R/3.4.1-gccmkl
module load R/3.4.1-gccmkl_20181025
sh ./sc_LinkData.sh D27_FACS
Rscript ../r.scripts/sc-TissueMapper_RUN.D27_FACS.R
bash.scripts/sc_TissueMapper-DS_D17.sh
+ 3
− 3
View file @ 142e5028
#!/bin/bash
#SBATCH --job-name R_FullAnalysis
#SBATCH -p 256GB,256GBv1
#SBATCH --job-name R_FullAnalysis.DS
#SBATCH -p 256GB,256GBv1,384GB
#SBATCH -N 1
#SBATCH -t 7-0:0:0
#SBATCH -o job_%j.out
......@@ -8,6 +8,6 @@
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load R/3.4.1-gccmkl
module load R/3.4.1-gccmkl_20181025
Rscript ../r.scripts/sc-TissueMapper_RUN.DS_D17.R
bash.scripts/sc_TissueMapper-Pd.sh
+ 6
− 3
View file @ 142e5028
#!/bin/bash
#SBATCH --job-name R_FullAnalysis
#SBATCH --job-name R_FullAnalysis.Pd
#SBATCH -p 256GB,256GBv1,384GB
#SBATCH -N 1
#SBATCH -t 7-0:0:0
......@@ -8,6 +8,9 @@
#SBATCH --mail-type ALL
#SBATCH --mail-user gervaise.henry@utsouthwestern.edu
module load R/3.4.1-gccmkl
module load R/3.4.1-gccmkl_20181025
Rscript ../r.scripts/sc-TissueMapper_RUN.Pd.R
\ No newline at end of file
sh ./sc_LinkData.sh Pd
Rscript ../r.scripts/sc-TissueMapper_RUN.Pd.R
Rscript ../r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R
genesets/genes.deg.Epi.csv
+ 769
− 43
View file @ 142e5028
This diff is collapsed.
genesets/genes.deg.LE.csv
+ 71
− 79
View file @ 142e5028
"","p_val","avg_logFC","pct.1","pct.2","p_val_adj"
"MSMB",0,4.37447239718553,0.686,0.25,0
"KLK3",0,3.78840323642459,0.672,0.08,0
"ACPP",0,3.13382648140614,0.596,0.031,0
"PLA2G2A",0,2.81892272901388,0.444,0.02,0
"KLK2",0,2.49353927858673,0.594,0.025,0
"TFF1",0,2.15334059243419,0.329,0.041,0
"AZGP1",0,2.0049490413152,0.897,0.345,0
"MT1G",0,1.98581083516144,0.254,0.053,0
"TSPAN8",0,1.82868329126073,0.741,0.077,0
"CLDN3",0,1.81872600017444,0.678,0.078,0
"NEAT1",0,1.66342555320245,0.992,0.932,0
"SEC11C",0,1.50139213870143,0.616,0.159,0
"RDH11",0,1.43938695807039,0.54,0.1,0
"CPE",0,1.42726881338742,0.609,0.195,0
"AGR2",0,1.38739010258592,0.569,0.185,0
"ADIRF",0,1.20378597343587,0.662,0.266,0
"TMPRSS2",0,1.17657374095478,0.598,0.1,0
"CD24",0,1.169691225058,0.679,0.256,0
"NPDC1",0,1.09945973647627,0.691,0.197,0
"ZG16B",0,1.09632759033173,0.705,0.381,0
"DHRS7",0,1.0780600230455,0.538,0.175,0
"SLC12A2",0,1.05453438337529,0.48,0.056,0
"ALDH1A3",0,1.00938272421439,0.708,0.409,0
"SORD",0,1.00762400392571,0.38,0.101,0
"KLK4",0,1.00075010108915,0.383,0.007,0
"SMIM14",0,0.990064601095013,0.664,0.201,0
"LMO7",0,0.981449368568302,0.611,0.256,0
"TMC5",0,0.944131304887237,0.446,0.062,0
"HMGCS1",0,0.928927221396835,0.375,0.112,0
"TPD52",0,0.921812891053796,0.65,0.253,0
"PPP3CA",0,0.918524421978158,0.497,0.095,0
"RASEF",0,0.918474385055297,0.443,0.041,0
"SLC44A4",0,0.915536056294867,0.565,0.098,0
"CIB1",0,0.910802697071163,0.822,0.534,0
"DSC2",0,0.908460413879168,0.476,0.123,0
"ELF3",0,0.90825329652683,0.83,0.543,0
"PART1",0,0.897849167503269,0.422,0.041,0
"PSCA",0,0.883834790597298,0.406,0.073,0
"INSIG1",0,0.879088049701038,0.435,0.138,0
"UQCRQ",0,0.851514992774378,0.826,0.57,0
"SLC39A6",0,0.847303216462186,0.565,0.186,0
"SERF2",0,0.837725135431806,0.997,0.944,0
"B2M",0,0.829136712088592,0.999,0.992,0
"SPINK1",0,0.823833574447547,0.288,0.064,0
"FXYD3",0,0.795341516329355,0.919,0.751,0
"STEAP2",0,0.786976946361527,0.386,0.061,0
"ABHD2",0,0.779715640605612,0.533,0.159,0
"CDC42EP5",0,0.774674582069486,0.613,0.217,0
"RAB27B",0,0.766534607182507,0.37,0.038,0
"GOLM1",0,0.760306018074098,0.393,0.037,0
"LINC00844",0,0.755952806149427,0.266,0.024,0
"LIPH",0,0.755836349227058,0.363,0.06,0
"MT1F",0,0.735175225812285,0.275,0.06,0
"BACE2",0,0.734623852847882,0.652,0.299,0
"FLJ20021",0,0.722227423027725,0.301,0.047,0
"ELOVL5",0,0.713313331698034,0.567,0.226,0
"LURAP1L",0,0.701460739959854,0.38,0.098,0
"ERBB3",0,0.701087186907271,0.429,0.074,0
"RNF144B",1.55475455202509e-306,0.760182455666201,0.38,0.119,3.61480433345833e-302
"GDF15",2.01211114811691e-304,1.02000284299098,0.482,0.18,4.67815841937181e-300
"MARCKSL1",3.87076430258232e-298,0.770455303625736,0.608,0.273,8.99952700350389e-294
"DBI",1.31508887521865e-297,0.859804208266948,0.708,0.443,3.05758163488335e-293
"ANKRD36C",1.31868558656178e-293,0.835820434758738,0.393,0.13,3.06594398875615e-289
"IL6ST",3.47797112066833e-286,0.72159960498644,0.532,0.237,8.08628285555386e-282
"TSC22D1",8.48846185126225e-286,0.981143608106769,0.72,0.435,1.97356738041847e-281
"H2AFJ",1.02086674628759e-266,0.715945859194024,0.799,0.556,2.37351518511864e-262
"NDRG1",1.60403620345049e-258,0.831848305281009,0.566,0.291,3.72938417302238e-254
"PLPP1",3.7845405697958e-257,0.852158040389265,0.364,0.126,8.79905682477522e-253
"SAT1",1.77156969700026e-217,0.792072293664311,0.988,0.968,4.1188995455256e-213
"TXNIP",1.17916945396995e-214,0.824032953148113,0.559,0.3,2.74156898048013e-210
"RPS27L",3.3815667813779e-198,0.802242921285638,0.675,0.444,7.86214276670362e-194
"NKX3-1",4.7758881227118e-173,0.714910742253142,0.291,0.107,1.11039398853049e-168
"GLRX",1.82073919060387e-148,0.705291464806051,0.296,0.119,4.233218618154e-144
"MT1E",1.64237321428792e-145,1.85089781749558,0.454,0.266,3.81851772321941e-141
"MGST1",1.23115980748549e-121,0.741661437165127,0.481,0.29,2.86244655240376e-117
"FDPS",4.05857338821562e-111,0.779405255968596,0.412,0.229,9.43618312760133e-107
"GADD45B",2.47436197380962e-62,0.957716958200775,0.512,0.375,5.75289158910737e-58
"MT2A",4.48785503840556e-60,1.21424445779221,0.712,0.623,1.04342629642929e-55
"JUN",6.73977419884717e-60,0.807603248779742,0.684,0.58,1.56699750123197e-55
"MSMB",0,4.41464597666458,0.686,0.239,0
"KLK3",0,3.79927069647092,0.672,0.077,0
"ACPP",0,3.1386912452006,0.596,0.031,0
"PLA2G2A",0,2.82011634326759,0.444,0.02,0
"KLK2",0,2.49685405203231,0.594,0.024,0
"TFF1",0,2.14654338516735,0.329,0.049,0
"AZGP1",0,2.01859676603558,0.897,0.363,0
"MT1G",0,1.99384816623504,0.254,0.051,0
"TSPAN8",0,1.84829135427538,0.741,0.07,0
"CLDN3",0,1.52902244520835,0.678,0.129,0
"NEAT1",0,1.51900228702682,0.992,0.94,0
"SEC11C",0,1.50870296196156,0.616,0.162,0
"CPE",0,1.46883854734414,0.609,0.176,0
"RDH11",0,1.43791145325306,0.54,0.106,0
"AGR2",0,1.20104069415542,0.569,0.26,0
"TMPRSS2",0,1.11615430581965,0.598,0.134,0
"ZG16B",0,1.1068824310419,0.705,0.394,0
"NPDC1",0,1.10185293373806,0.691,0.205,0
"DHRS7",0,1.09703119088929,0.538,0.17,0
"SLC12A2",0,1.06008931420528,0.48,0.057,0
"CD24",0,1.0023510757434,0.679,0.287,0
"KLK4",0,1.0021842053216,0.383,0.007,0
"SMIM14",0,0.984221962278597,0.664,0.215,0
"ADIRF",0,0.946372464713426,0.662,0.323,0
"TMC5",0,0.935364519656032,0.446,0.069,0
"PPP3CA",0,0.922202053364421,0.497,0.099,0
"TPD52",0,0.913208727296692,0.65,0.27,0
"PART1",0,0.898217216909069,0.422,0.043,0
"DSC2",0,0.895541403843838,0.476,0.136,0
"CIB1",0,0.886862374771338,0.822,0.563,0
"SLC44A4",0,0.878684659118734,0.565,0.125,0
"RASEF",0,0.86881880871108,0.443,0.066,0
"INSIG1",0,0.85742625724725,0.435,0.154,0
"SLC39A6",0,0.857043457571851,0.565,0.189,0
"UQCRQ",0,0.834578123494899,0.826,0.593,0
"STEAP2",0,0.792982361190698,0.386,0.061,0
"SERF2",0,0.781333597608586,0.997,0.951,0
"RNF144B",0,0.773755237410753,0.38,0.117,0
"ABHD2",0,0.770356389121221,0.533,0.173,0
"LINC00844",0,0.763189426201448,0.266,0.022,0
"GOLM1",0,0.756286404275184,0.393,0.042,0
"RAB27B",0,0.753392949985878,0.37,0.048,0
"BACE2",0,0.744348788707172,0.652,0.307,0
"MT1F",0,0.738696732268015,0.275,0.061,0
"PSCA",0,0.736863698556448,0.406,0.126,0
"ELOVL5",0,0.735734416728113,0.567,0.224,0
"FLJ20021",0,0.717428972989397,0.301,0.052,0
"B2M",0,0.710373668966736,0.999,0.993,0
"PMEPA1",0,0.696773049825896,0.585,0.265,0
"HMGCS1",3.42461756531418e-303,0.921004561752512,0.375,0.123,7.96223583935546e-299
"LMO7",2.71147973468059e-301,0.864726032278975,0.611,0.306,6.30419038313236e-297
"SORD",8.31601295337064e-295,0.953794980653706,0.38,0.13,1.93347301165867e-290
"FXYD3",1.82555434125871e-291,0.72276053621229,0.919,0.777,4.24441384342649e-287
"ALDH1A3",2.81626138046117e-290,0.96151656515069,0.708,0.445,6.54780770957222e-286
"IL6ST",1.76482583457318e-279,0.720964598078865,0.532,0.248,4.10322006538265e-275
"TSC22D1",2.00482551149997e-278,0.973848388488826,0.72,0.451,4.66121931423742e-274
"DBI",4.58584954232928e-275,0.83482691237847,0.708,0.472,1.06621001859156e-270
"PLPP1",4.11040262169908e-270,0.863079043578098,0.364,0.126,9.55668609545036e-266
"H2AFJ",1.13241110585981e-266,0.716349174750414,0.799,0.574,2.63285582112405e-262
"ELF3",5.4285571136693e-254,0.750384596642815,0.83,0.589,1.26213952892811e-249
"SPINK1",4.00680293998379e-247,0.76488254323365,0.288,0.081,9.31581683546232e-243
"NDRG1",5.11762466241304e-211,0.773387403674286,0.566,0.327,1.18984773401103e-206
"RPS27L",2.73209808270906e-210,0.820414918744851,0.675,0.447,6.35212804229856e-206
"SAT1",4.68375548904372e-209,0.782698466152202,0.988,0.972,1.08897315120266e-204
"NKX3-1",1.05503565393114e-206,0.73914604464308,0.291,0.098,2.4529578953899e-202
"TXNIP",4.11598983981011e-200,0.81506575663786,0.559,0.317,9.56967637755851e-196
"MT1E",1.27120095406914e-171,1.90413450179974,0.454,0.251,2.95554221821074e-167
"FDPS",9.5545509588763e-102,0.772645903487662,0.412,0.241,2.22143309793874e-97
"MT2A",6.29978747466877e-90,1.3228614857529,0.712,0.597,1.46470058786049e-85
"JUN",1.13685849258666e-64,0.829884133857674,0.684,0.586,2.64319599526398e-60
"GADD45B",1.11647864383533e-59,0.960832459862578,0.512,0.387,2.59581284691715e-55
genesets/genes.deg.St.csv
+ 415
− 109
View file @ 142e5028
This diff is collapsed.
r.scripts/sc-TissueMapper_RUN.D17_FACS.R
+ 3
− 4
View file @ 142e5028
......@@ -167,13 +167,12 @@ if (opt$st==TRUE){
sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=0.5,folder="post.stress",red="pca")
}
gene.set1 <- read_delim("./genesets/genes.deg.Epi.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Epi"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.St.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set,gene.set1)
......
r.scripts/sc-TissueMapper_RUN.D27_FACS.R
+ 3
− 4
View file @ 142e5028
......@@ -167,13 +167,12 @@ if (opt$st==TRUE){
sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=0.5,folder="post.stress",red="pca")
}
gene.set1 <- read_delim("./genesets/genes.deg.Epi.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Epi"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.St.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set,gene.set1)
......
r.scripts/sc-TissueMapper_RUN.DS_D17.R
+ 46
− 6
View file @ 142e5028
......@@ -13,6 +13,8 @@ library(viridis)
library(reshape2)
library(NMI)
options(bitmapType="cairo")
source("../r.scripts/sc-TissueMapper.R")
#Create folder structure
......@@ -101,7 +103,7 @@ option_list=list(
make_option("--lx",action="store",default=0.2,type='numeric',help="x low threshold for hvg selection"),
make_option("--hx",action="store",default=5,type='numeric',help="x high threshold for hvg selection"),
make_option("--ly",action="store",default=1,type='numeric',help="y low threshold for hvg selection"),
make_option("--cc",action="store",default=TRUE,type='logical',help="Scale cell cycle?"),
make_option("--cc",action="store",default=FALSE,type='logical',help="Scale cell cycle?"),
make_option("--cca",action="store",default=50,type='integer',help="Number of CCAs to cacluate"),
make_option("--acca",action="store",default=30,type='integer',help="Number of CCAs to align"),
make_option("--pc",action="store",default=50,type='integer',help="Number of PCs to cacluate"),
......@@ -178,13 +180,11 @@ if (opt$st==TRUE){
sc10x <- scCluster(sc10x,pc.use=pc.use.poststress,res.use=opt$res.poststress,folder="post.stress",red="pca")
}
gene.set1 <- read_delim("./genesets/genes.deg.Epi.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Epi"
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.St.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "St"
gene.set <- c(gene.set,gene.set1)
......@@ -305,6 +305,11 @@ rm(gene.set)
sc10x.Epi.NE <- scNE(sc10x.Epi,neg="dws",cut=opt$cut.ne)
sc10x <- scMerge(sc10x,sc10x.Epi,sc10x.St,i.1="Epi.dws.sc",i.2="St.dws.sc",nm="Merge_Epi.dws.sc_St.dws.sc")
cells.ne <- names(sc10x.Epi.NE@ident[sc10x.Epi.NE@ident=="NE"])
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws.sc_St.dws.sc")
sc10x <- SetIdent(object=sc10x,cells.use=cells.ne,ident.use="NE")
sc10x <- StashIdent(object=sc10x,save.name="Merge_Epi.dws.sc_St.dws.sc_NE")
rm(cells.ne)
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws.sc_St.dws.sc")
sc10x <- SetAllIdent(object=sc10x,id="Merge_Epi.dws.sc_St.dws.sc")
......@@ -316,7 +321,7 @@ plot <- plot+guides(colour=guide_legend(override.aes=list(size=10)))
plot(plot)
dev.off()
scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc")
scTables(sc10x,i.1="samples",i.2="Merge_Epi.dws.sc_St.dws.sc_NE")
sctSNECustCol(sc10x,i="Lin",bl="Epi",rd="St",file="D17")
sctSNECustCol(sc10x,i="Merge_Epi.dws.sc_St.dws.sc",bl=c("BE","LE","OE_SCGB","OE_KRT13"),rd=c("Fib","SM","Endo","Leu"),file="D17")
......@@ -344,6 +349,41 @@ for (g in c("Fib","SM","Endo","Leu")){
rm(i)
rm(g)
try(genes.deg.Stress <- scDEG(sc10x.Stress,i="Stress",g.1="Stress",g.2="ALL",pct=0.5,t=5))
try(genes.deg.Epi <- scDEG(sc10x,i="Lin",g.1="Epi",g.2="St",t=2))
try(genes.deg.St <- scDEG(sc10x,i="Lin",g.1="St",g.2="Epi",t=2))
try(genes.deg.BE <- scDEG(sc10x,i="Merge_Epi.dws.sc_St.dws.sc_NE",g.1="BE",g.2=c("LE","OE_SCGB","OE_KRT13"),pct=0.25,t=2))
try(genes.deg.LE <- scDEG(sc10x,i="Merge_Epi.dws.sc_St.dws.sc_NE",g.1="LE",g.2=c("BE","OE_SCGB","OE_KRT13"),pct=0.25,t=2))
try(genes.deg.OE1 <- scDEG(sc10x,i="Merge_Epi.dws.sc_St.dws.sc_NE",g.1="OE_SCGB",g.2=c("BE","LE","OE_KRT13"),pct=0.25,t=2))
try(genes.deg.OE2 <- scDEG(sc10x,i="Merge_Epi.dws.sc_St.dws.sc_NE",g.1="OE_KRT13",g.2=c("BE","LE","OE_SCGB"),pct=0.25,t=2))
try(genes.deg.NE <- scDEG(sc10x,i="Merge_Epi.dws.sc_St.dws.sc_NE",g.1="NE",g.2=c("BE","LE","OE_SCGB","OE_KRT13"),pct=0.01,t=1))
try(genes.deg.Fib <- scDEG(sc10x.St,i="Merge_Epi.dws.sc_St.dws.sc",g.1="Fib",g.2=c("SM","Endo","Leu"),pct=0.25,t=2))
try(genes.deg.SM <- scDEG(sc10x.St,i="Merge_Epi.dws.sc_St.dws.sc",g.1="SM",g.2=c("Fib","Endo","Leu"),pct=0.25,t=2))
try(genes.deg.Endo <- scDEG(sc10x.St,i="Merge_Epi.dws.sc_St.dws.sc",g.1="Endo",g.2=c("Fib","SM","Leu"),pct=0.25,t=2))
try(genes.deg.Leu <- scDEG(sc10x.St,i="St.dws.sc",g.1="Leu",g.2=c("Fib","SM","Endo"),pct=0.25,t=2))
try(genes.deg.BE.unique <- setdiff(rownames(genes.deg.BE),Reduce(union,list(rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.LE.unique <- setdiff(rownames(genes.deg.LE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.OE1.unique <- setdiff(rownames(genes.deg.OE1),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.OE2.unique <- setdiff(rownames(genes.deg.OE2),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.NE.unique <- setdiff(rownames(genes.deg.NE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.Fib.unique <- setdiff(rownames(genes.deg.Fib),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.SM.unique <- setdiff(rownames(genes.deg.SM),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.Endo),rownames(genes.deg.Leu)))))
try(genes.deg.Endo.unique <- setdiff(rownames(genes.deg.Endo),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Leu)))))
try(genes.deg.Leu.unique <- setdiff(rownames(genes.deg.Leu),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo)))))
try(genes.deg.5 <- c(genes.deg.BE.unique[1:5],genes.deg.LE.unique[1:5],genes.deg.OE1.unique[1:5],genes.deg.OE2.unique[1:5],genes.deg.NE.unique[1:5],genes.deg.Fib.unique[1:5],genes.deg.SM.unique[1:5],genes.deg.Endo.unique[1:5],genes.deg.Leu.unique[1:5]))
try(genes.deg.5 <- rev(genes.deg.5))
try(genes.deg.10 <- c(genes.deg.BE.unique[1:10],genes.deg.LE.unique[1:10],genes.deg.OE1.unique[1:10],genes.deg.OE2.unique[1:10],genes.deg.NE.unique[1:10],genes.deg.Fib.unique[1:10],genes.deg.SM.unique[1:10],genes.deg.Endo.unique[1:10],genes.deg.Leu.unique[1:10]))
try(genes.deg.10 <- rev(genes.deg.10))
for (i in ls(pattern="^genes.deg")){
try(write.table(get(i),file=paste0("./analysis/deg/",i,".csv"),row.names=TRUE,col.names=NA,append=FALSE,sep=","))
}
save(list=ls(pattern="sc10x.Stress"),file="./analysis/sc10x.Stress.Rda")
rm(list=ls(pattern="sc10x.Stress"))
......
r.scripts/sc-TissueMapper_RUN.DS_D17.aggr.R
+ 2
− 0
View file @ 142e5028
......@@ -13,6 +13,8 @@ library(viridis)
library(reshape2)
library(NMI)
options(bitmapType="cairo")
source("../r.scripts/sc-TissueMapper.R")
setwd("../")
......
r.scripts/sc-TissueMapper_RUN.Pd.R
+ 9
− 3
View file @ 142e5028
......@@ -11,6 +11,8 @@ library(monocle)
library(dplyr)
library(viridis)
options(bitmapType="cairo")
source("../r.scripts/sc-TissueMapper.R")
#Create folder structure
......@@ -407,14 +409,16 @@ scTables(sc10x,i.1="Merge_Epi.dws_St.go_NE",i.2="Merge_Epi.dws_St.go")
genes.deg.Stress <- scDEG(sc10x.Stress,i="Stress",g.1="Stress",g.2="ALL",pct=0.5,t=5)
genes.deg.Epi <- scDEG(sc10x,i="Lin",g.1="Epi",g.2="St",t=2)
genes.deg.St <- scDEG(sc10x,i="Lin",g.1="St",g.2="Epi",t=2)
genes.deg.Epi <- scDEG(sc10x,i="Lin",g.1="Epi",g.2="St",t=0)
genes.deg.St <- scDEG(sc10x,i="Lin",g.1="St",g.2="Epi",t=0)
genes.deg.BE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="BE",g.2=c("LE","OE1","OE2"),pct=0.25,t=2)
genes.deg.LE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="LE",g.2=c("BE","LE","OE1"),pct=0.25,t=2)
genes.deg.LE <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="LE",g.2=c("BE","OE1","OE2"),pct=0.25,t=2)
genes.deg.OE1 <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="OE1",g.2=c("BE","LE","OE2"),pct=0.25,t=2)
genes.deg.OE2 <- scDEG(sc10x.Epi.NE,i="Epi.dws.sub",g.1="OE2",g.2=c("BE","LE","OE1"),pct=0.25,t=2)
genes.deg.Epi <- genes.deg.Epi[rownames(genes.deg.Epi) %in% setdiff(rownames(genes.deg.Epi),union(union(union(rownames(genes.deg.BE),rownames(genes.deg.LE)),rownames(genes.deg.OE1)),rownames(genes.deg.OE2))),]
genes.deg.NE <- scDEG(sc10x.Epi.NE,i="NE",g.1="NE",g.2="ALL",pct=0.01,t=1)
genes.deg.Fib <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="Fib",g.2=c("SM","Endo","Leu"),pct=0.25,t=2)
......@@ -422,6 +426,8 @@ genes.deg.SM <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="SM",g.2=c("Fib","
genes.deg.Endo <- scDEG(sc10x.St,i="Merge_Epi.dws_St.go_NE",g.1="Endo",g.2=c("Fib","SM","Leu"),pct=0.25,t=2)
genes.deg.Leu <- scDEG(sc10x.St,i="St.go",g.1="Leu",g.2=c("Fib","SM","Endo"),pct=0.25,t=2)
genes.deg.St <- genes.deg.St[rownames(genes.deg.St) %in% setdiff(rownames(genes.deg.St),union(union(union(rownames(genes.deg.Fib),rownames(genes.deg.SM)),rownames(genes.deg.Endo)),rownames(genes.deg.Leu))),]
genes.deg.BE.unique <- setdiff(rownames(genes.deg.BE),Reduce(union,list(rownames(genes.deg.LE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.LE.unique <- setdiff(rownames(genes.deg.LE),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.OE1),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
genes.deg.OE1.unique <- setdiff(rownames(genes.deg.OE1),Reduce(union,list(rownames(genes.deg.BE),rownames(genes.deg.LE),rownames(genes.deg.OE2),rownames(genes.deg.NE),rownames(genes.deg.Fib),rownames(genes.deg.SM),rownames(genes.deg.Endo),rownames(genes.deg.Leu))))
......
r.scripts/sc-TissueMapper_RUN.Pd.StressCompare.R 0 → 100755
+ 135
− 0
View file @ 142e5028
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(VennDiagram)
library(Cairo)
options(bitmapType="cairo")
if (!dir.exists("../analysis/stress")){
dir.create("../analysis/stress")
}
load("../analysis/sc10x.Stress.Rda")
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL")
gene.exp <- as.data.frame(as.matrix(sc10x.Stress@raw.data[,colnames(sc10x.Stress@data)]))
genes.stress.go <- read_csv("../genesets/DEG_C2.CGP.M10970.txt")
genes.stress.go <- genes.stress.go[2:nrow(genes.stress.go),]
genes.stress.go <- as.list(genes.stress.go)
genes.stress.dws <- read_delim("../genesets/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress.dws <- genes.stress.dws[1]
colnames(genes.stress.dws) <- "scDWS.Stress"
genes.stress.dws <- as.list(genes.stress.dws)
#genes.stress.vdb <- read_delim("./vanderBrink.Stress.mus.tsv","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
#colnames(genes.stress.vdb) <- "vdBrink.Stress"
genes.stress.vdb <- c("Actg1__chr11","Ankrd1__chr19","Arid5a__chr1","Atf3__chr1","Atf4__chr15","Bag3__chr7","Bhlhe40__chr6","Brd2__chr17","Btg1__chr10","Btg2__chr1","Ccnl1__chr3","Ccrn4l__chr3","Cebpb__chr2","Cebpd__chr16","Cebpg__chr7","Csrnp1__chr9","Cxcl1__chr5","Cyr61__chr3","Dcn__chr10","Ddx3x__chrX","Ddx5__chr11","Des__chr1","Dnaja1__chr4","Dnajb1__chr8","Dnajb4__chr3","Dusp1__chr17","Dusp8__chr7","Egr1__chr18","Egr2__chr10","Eif1__chr11","Eif5__chr12","Erf__chr7","Errfi1__chr4","Fam132b__chr1","Fos__chr12","Fosb__chr7","Fosl2__chr5","Gadd45a__chr6","Gcc1__chr6","Gem__chr4","H3f3b__chr11","Hipk3__chr2","Hsp90aa1__chr12","Hsp90ab1__chr17","Hspa1a__chr17","Hspa1b__chr17","Hspa5__chr2","Hspa8__chr9","Hspb1__chr5","Hsph1__chr5","Id3__chr4","Idi1__chr13","Ier2__chr8","Ier3__chr17","Ifrd1__chr12","Il6__chr5","Irf1__chr11","Irf8__chr8","Itpkc__chr7","Jun__chr4","Junb__chr8","Jund__chr8","Klf2__chr8","Klf4__chr4","Klf6__chr13","Klf9__chr19","Litaf__chr16","Lmna__chr3","Maff__chr15","Mafk__chr5","Mcl1__chr3","Midn__chr10","Mir22hg__chr11","Mt1__chr8","Mt2__chr8","Myadm__chr7","Myc__chr15","Myd88__chr9","Nckap5l__chr15","Ncoa7__chr10","Nfkbia__chr12","Nfkbiz__chr16","Nop58__chr1","Nppc__chr1","Nr4a1__chr15","Odc1__chr12","Osgin1__chr8","Oxnad1__chr14","Pcf11__chr7","Pde4b__chr4","Per1__chr11","Phlda1__chr10","Pnp__chr14","Pnrc1__chr4","Ppp1cc__chr5","Ppp1r15a__chr7","Pxdc1__chr13","Rap1b__chr10","Rassf1__chr9","Rhob__chr12","Rhoh__chr5","Ripk1__chr13","Sat1__chrX","Sbno2__chr10","Sdc4__chr2","Serpine1__chr5","Skil__chr3","Slc10a6__chr5","Slc38a2__chr15","Slc41a1__chr1","Socs3__chr11","Sqstm1__chr11","Srf__chr17","Srsf5__chr12","Srsf7__chr17","Stat3__chr11","Tagln2__chr1","Tiparp__chr3","Tnfaip3__chr10","Tnfaip6__chr2","Tpm3__chr3","Tppp3__chr8","Tra2a__chr6","Tra2b__chr16","Trib1__chr15","Tubb4b__chr2","Tubb6__chr18","Ubc__chr5","Usp2__chr9","Wac__chr18","Zc3h12a__chr4","Zfand5__chr19","Zfp36__chr7","Zfp36l1__chr12","Zfp36l2__chr17","Zyx__chr6","Gadd45g__chr13","Hspe1__chr1","Ier5__chr1","Kcne4__chr1")
genes.stress.vdb <- gsub("\\__.*","",genes.stress.vdb)
genes.stress.vdb <- data.frame(vdBrink.Stress=genes.stress.vdb)
genes.homolog <- read_delim("../genesets/Ensemble.mus-hum.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress.vdb <- list("vdBrink.Stress"=merge(genes.stress.vdb,genes.homolog[,c(2,4)],by.x="vdBrink.Stress",by.y="Gene name")[,2])
venn.diagram(list(go=unlist(genes.stress.go),dws=unlist(genes.stress.dws),vdb=unlist(genes.stress.vdb)),"../analysis/stress/Venn.tiff")
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
gene.stress.pct <- apply(gene.exp,2,function(x) sum(x[rownames(gene.exp) %in% unlist(genes.stress)])/sum(x))
#histo <- hist(gene.stress.pct,breaks=100,prob=TRUE,plot=TRUE,main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes")
#abline(v=.0575,col="red")
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL+Stress")
cell.index.stress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="Stress"])
cell.index.stress <- which(colnames(gene.exp) %in% cell.index.stress)
cell.index.ntstress <- names(sc10x.Stress@ident[sc10x.Stress@ident=="ALL"])
cell.index.ntstress <- which(colnames(gene.exp) %in% cell.index.ntstress)
histo.stress <- hist(gene.stress.pct[cell.index.stress],breaks=seq(0,ceiling(max(gene.stress.pct)*100)/100,0.005),main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes",plot=FALSE)
#abline(v=.0575,col="red")
histo.ntstress <- hist(gene.stress.pct[cell.index.ntstress],breaks=seq(0,ceiling(max(gene.stress.pct)*100)/100,0.005),main="Distribution of Stress Gene Expression in Transcriptome",xlab="% of transcriptome associated to stress genes",plot=FALSE)
#abline(v=.0575,col="red")
Cairo(paste0("../analysis/stress/Hist.Stress.",i,".png"),width=1000,height=500,bg="white")
plot(histo.ntstress$mids,histo.ntstress$density,col=rgb(0,1,0,1/5),xlim=c(0,ceiling(max(gene.stress.pct)*100)/100),ylim=c(0,ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10),type="h",lend="square",lwd=10,xlab="Fraction of transcriptome associated to stress genes",ylab="Density")
par(new=TRUE)
plot(histo.stress$mids,histo.stress$density,col=rgb(1,0,0,1/5),xlim=c(0,ceiling(max(gene.stress.pct)*100)/100),ylim=c(0,ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10),type="h",lend="square",lwd=10,axes=FALSE,xlab="",ylab="")
abline(v=0.0575,col="red")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks>0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks>0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,paste0(round(sum(histo.ntstress$counts[histo.ntstress$breaks<=0.035][!is.na(histo.ntstress$counts[histo.ntstress$breaks<=0.035])])/sum(histo.ntstress$counts[!is.na(histo.ntstress$counts)])*100),"%"),col="green")
text(0.0575+0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks>0.035][!is.na(histo.stress$counts[histo.stress$breaks>0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575-0.01,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,paste0(round(sum(histo.stress$counts[histo.stress$breaks<=0.035][!is.na(histo.stress$counts[histo.stress$breaks<=0.035])])/sum(histo.stress$counts[!is.na(histo.stress$counts)])*100),"%"),col="red")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3+2.5,"Not Stressed",col="green")
text(0.0575+0.05,(ceiling(max(histo.stress$density,histo.ntstress$density)/10)*10)/4*3-2.5,"Stressed",col="red")
dev.off()
}
rm(i)
rm(genes.stress)
rm(gene.stress.pct)
rm(cell.index.stress)
rm(cell.index.ntstress)
rm(histo.stress)
rm(histo.ntstress)
for (cut in c(0.9,0.95)){
for (i in c("go","dws","vdb")){
genes.stress=get(paste0("genes.stress.",i))
sc10x.Stress <- SetAllIdent(object=sc10x.Stress,id="ALL")
sc10x.Stress.dta <- as.data.frame(as.matrix(sc10x.Stress@data[rownames(sc10x.Stress@scale.data) %in% unlist(genes.stress),]))
#Run PCA of subsetted stress genes
sc10x.Stress.dta <- t(sc10x.Stress.dta)
sc10x.Stress.dta <- sc10x.Stress.dta[,apply(sc10x.Stress.dta,2,var)!=0]
sc10x.Stress.dta.pca <- prcomp(sc10x.Stress.dta,center=TRUE,scale.=TRUE)
sc10x.Stress.dta.pca <- sc10x.Stress.dta.pca$x[,1:2]
colnames(x=sc10x.Stress.dta.pca) <- paste0("Stress",1:2)
if (skewness(sc10x.Stress.dta.pca[,1])<0){
sc10x.Stress.dta.pca[,1] <- (-sc10x.Stress.dta.pca[,1])
}
if (skewness(sc10x.Stress.dta.pca[,2])<0){
sc10x.Stress.dta.pca[,2] <- (-sc10x.Stress.dta.pca[,2])
}
sc10x.dta <- sc10x.Stress
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="cell.embeddings",new.data=sc10x.Stress.dta.pca)
sc10x.dta <- SetDimReduction(object=sc10x.dta,reduction.type="Stress",slot="key",new.data="Stress")
#CDS
cdf <- ecdf(GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
cut.x <- quantile(cdf,probs=cut)
assign(i,GetCellEmbeddings(object=sc10x.dta,reduction.type="Stress",dims.use=1))
assign(paste0("cut.",i,".",cut),cut.x)
}}
rm(cut)
rm(i)
rm(genes.stress)
rm(sc10x.Stress.dta)
rm(sc10x.Stress.dta.pca)
rm(sc10x.dta)
rm(cdf)
rm(cut.x)
comp <- cbind(vdb,go[match(rownames(vdb),rownames(go))],dws[match(rownames(vdb),rownames(dws))])
colnames(comp) <- c("vdb","go","dws")
comp <- data.frame(comp)
Cairo(paste0("../analysis/stress/Stress.vdb.vs.go.png"),width=1000,height=500,bg="white")
ggplot(comp,aes(x=vdb,y=go))+geom_point(size=0.01)+geom_vline(xintercept=cut.vdb.0.95)+geom_hline(yintercept=cut.go.0.95)+annotate("text",x=max(comp$vdb)/4*3,y=max(comp$go)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb>cut.vdb.0.95,][comp[comp$vdb>cut.vdb.0.95,"go"]>cut.go.0.95,])/nrow(comp[comp$vdb>cut.vdb.0.95,])*100,1),"%"))+annotate("text",x=max(comp$vdb)/100,y=max(comp$go)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb<=cut.vdb.0.95,][comp[comp$vdb<=cut.vdb.0.95,"go"]<=cut.go.0.95,])/nrow(comp[comp$vdb<=cut.vdb.0.95,])*100,1),"%"))
dev.off()
Cairo(paste0("../analysis/stress/Stress.vdb.vs.dws.png"),width=1000,height=500,bg="white")
ggplot(comp,aes(x=vdb,y=dws))+geom_point(size=0.01)+geom_vline(xintercept=cut.vdb.0.95)+geom_hline(yintercept=cut.dws.0.95)+annotate("text",x=max(comp$vdb)/4*3,y=max(comp$dws)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb>cut.vdb.0.95,][comp[comp$vdb>cut.vdb.0.95,"dws"]>cut.dws.0.95,])/nrow(comp[comp$vdb>cut.vdb.0.95,])*100,1),"%"))+annotate("text",x=max(comp$vdb)/100,y=max(comp$dws)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$vdb<=cut.vdb.0.95,][comp[comp$vdb<=cut.vdb.0.95,"dws"]<=cut.dws.0.95,])/nrow(comp[comp$vdb<=cut.vdb.0.95,])*100,1),"%"))
dev.off()
Cairo(paste0("../analysis/stress/Stress.go.vs.dws.png"),width=1000,height=500,bg="white")
ggplot(comp,aes(x=go,y=dws))+geom_point(size=0.01)+geom_vline(xintercept=cut.go.0.95)+geom_hline(yintercept=cut.dws.0.95)+annotate("text",x=max(comp$go)/4*3,y=max(comp$dws)/4*3,col="red",size=10,fontface=2,label=paste0(round(nrow(comp[comp$go>cut.go.0.95,][comp[comp$go>cut.go.0.95,"dws"]>cut.dws.0.95,])/nrow(comp[comp$go>cut.go.0.95,])*100,1),"%"))+annotate("text",x=max(comp$go)/100,y=max(comp$dws)/100,col="green",size=10,fontface=2,label=paste0(round(nrow(comp[comp$go<=cut.go.0.95,][comp[comp$go<=cut.go.0.95,"dws"]<=cut.dws.0.95,])/nrow(comp[comp$go<=cut.go.0.95,])*100,1),"%"))
dev.off()
rm(comp)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment