search

lanagarmire / Asgard

Other
37 stars 15 forks source link

DrugScore function return NAN and 0 #7

Open saisaitian opened 1 year ago

saisaitian commented 1 year ago

Hi,

Thanks for all the work on this tool. I am having trouble getting results for DrugScore function. I tried a few other tissue types including breast, kidney, lymphatic tissue and they seemed also return 0 and NAN in DrugScore function. image

My code is below

Drug repurposing using drugs/compounds

Gene.list<-readRDS("TNBC_genelist_limma.rds") my_gene_info<-read.table(file="DrugReference/skin_gene_info.txt",sep="\t",header = T,quote = "") my_drug_info<-read.table(file="DrugReference/skin_drug_info.txt",sep="\t",header = T,quote = "") cmap.ref.profiles<-GetDrugRef(drug.response.path = 'DrugReference/skin_rankMatrix.txt', probe.to.genes = my_gene_info, drug.info = my_drug_info) Drug.ident.res = GetDrug(gene.data = Gene.list, drug.ref.profiles = cmap.ref.profiles, repurposing.unit = "drug", connectivity = "negative", drug.type="all") saveRDS(Drug.ident.res,file="TNBC_drugs_limma_all.rds")

Drug/compound score

library(cmapR) SC.integrated<-readRDS("TNBC_SCdata.rds") Gene.list<-readRDS("TNBC_genelist_limma.rds") Drug.ident.res<-readRDS("TNBC_drugs_limma_all.rds") GSE92742.gctx.path="./Lincs/GSE92742_Broad_LINCS_Level5_COMPZ.MODZ_n473647x12328.gctx" GSE70138.gctx.path="./Lincs/GSE70138_Broad_LINCS_Level5_COMPZ_n118050x12328_2017-03-06.gctx" Tissue="breast"

Drug.score<-DrugScore(SC.integrated=SC.integrated, Gene.data=Gene.list, Cell.type=NULL, Drug.data=Drug.ident.res, FDA.drug.only=F, Case=Case, Tissue="breast", GSE92742.gctx=GSE92742.gctx.path, GSE70138.gctx=GSE70138.gctx.path)

hebinghb commented 1 year ago

I noticed you used the parameter Case, which defined the samples used for calculating drug score, but didn't see where you define it in your code. I'm not sure if it is the issue. But if you used the wrong Case names that do not exist, it may result in NAN or 0 output of drug score. Here is the Case defined in the demo code:

Case sample names

Case=c("PDX-110","PDX-332")

NarainritKaruna commented 9 months ago

Hi,

Tissue="muscle" Cell.type = c("B-cells","Cardiomyocytes","Endocardial cells", "Endothelial cells", "Fibroblasts", "Granulocytes", "Macrophages", "Pericytes", "Schwann cells", "Smooth muscle cells", "T-cells")

DbCM

Case <- c("HFD_1", "HFD_2", "HFD_3") Drug.score.bcell<-DrugScore(SC.integrated=SC.data, Gene.data=Gene.list, Cell.type= Cell.type, Drug.data=Drug.ident.res, FDA.drug.only=T, Case=Case, Tissue=Tissue, GSE92742.gctx=GSE92742.gctx.path, GSE70138.gctx=GSE70138.gctx.path)

I got the same problem for running all cell types together. image

Then, try to look into each cell type; the score is shown. I think it might the process of calculation when putting all cell types together.

Drug.score.mac<-DrugScore(SC.integrated=SC.data, Gene.data=Gene.list, Cell.type="Macrophages", Drug.data=Drug.ident.res, FDA.drug.only=T, Case=Case, Tissue=Tissue, GSE92742.gctx=GSE92742.gctx.path, GSE70138.gctx=GSE70138.gctx.path)

image

NarainritKaruna commented 9 months ago

I figure out how the function calculated NaN (usually from 0/0)

I found that when

Then, It will be 0/0 by Ratio.treated <- length(D.genes.treated)/length(D.D.genes)

for (Drug in C.Drugs) { D.genes.treated <- NULL drug.treatments <- subset(data_infor,pert_iname == Drug)$sig_id if(length(drug.treatments)>1){ drug.responses <- data[,drug.treatments] drug.responses.mean <- apply(drug.responses,1,mean) }else{ drug.responses <- data[,drug.treatments] drug.responses.mean <- drug.responses } D.D.genes <- intersect(names(D.gene.expression),names(drug.responses.mean)) D.genes.treated <- D.gene.expression[D.D.genes]drug.responses.mean[D.D.genes] D.genes.treated <- D.genes.treated[which(D.genes.treated>0)] Mean.treated <- mean(D.genes.treated) Ratio.treated <- length(D.genes.treated)/length(D.D.genes) Coverage.treated <- Drug.coverage[Drug]/100 Treated.score <- (Ratio.treatedCoverage.treated) Single.treated.score.list <- c(Single.treated.score.list,Treated.score)

Best,

hebinghb commented 9 months ago

Thanks for raising this issue. If drug.responses.mean got NULL, it came out because there is no treatment response data for the selected drugs in the selected L1000 tissue. If the both D.gene.expression and drug.responses.mean have value, but D.D.genes got NULL, it means there is no overlap between the disease-related differential genes and drug response genes.