This is part 2 to the previous quick analysis and predictive analytics testing in finding the top genes with fold change values from a study on nasopharyngeal carcinoma (NPC) and how Epstein-Barr viral infection (EBV) influences the proliferation rate to increase and to create more mutant strands of the target protein with an amino acid sequence change from Histadine to Arginine at the 101th amino acid marker of the messenger RNA or mRNA strand. I reviewed the pdf free study for EBVaNPC GOF and found that it is their term for Epstein-Barr Virus associated Gene Ontology Fingerprint. They found the GOF by screening thousands of studies on EBV and tested the differentially expressed genes then compared those to a network of gene sets from gene set expression analysis or GSEA for common genes as influencers to the type of NPC that is more aggressive or has seen a marked change from commercial line NPC they call the Wild Type (WT) vs their mutatated (MUT) strand that they used strand profilers to create the mutation at H101R that they found in their own clinical study of 26 EBV infected patients where 13 had EBV and the other 13 had EBV and NPC. They saw a pattern of changes in the Latent Membrane Protein (LMP1) to have the histadine amino acid (H amino acid letter) replacement at the 101th or 101st amino acid in the strand with Arginine (R amino acid letter).
Additonally in my summary of the report that was free to download at Pubmed article. The study did much lab and computational analysis and developed its own strategy for finding various lengths of strands that are created to use in single cell RNA sequencing so that they ranked the strands and found the length of the strand of nucleotides had negligible affects on the outcome of their own predictive machine learners to predict the class of their samples as EBV or EBV with NPC. The LMP1 gene is said to be an intracellular signaling protein and affects the T cell activity in identifying antigens that need removal or cell apoptosis of an unhealthy or infected cell. When there is a defective LMP1 gene product it causes an increase in a gene, FOXP3, that causes T cell anergy, looking it up this is the same process that makes HIV exhaust the T cells from using their CD4 proliferation to get natural killer cells in the bone marrow to (MEN - all men are natural killers - mnemonic for the Monocytes, Eosinophils, and Neutrophils that are made in bone marrow and only the Eosinophils and Neutrophils are granolocytes that send out granules to destroy infected cells due to other mnemonic BEN where Basophiles, Eosinophils, and Neutrophils are only White Blood Cells that have granulocytes, and monocytes are macrophages once they pass the blood vessels into the intracellular spaces of tissue in a response to injury defense producing 5 cardinal signs of inflammation or redness (rubor), pain (dolor), swelling (tumor), heat (calor), and funcio de laisso or temporary loss of function).
The mutated strand of H101R of the gene LMP1 in EBV infected patients can lead to NPCby increasing the FOXP3 gene that creates T-cell anergy or stops the immune regulating defenses of natural killer cells from seeking out and destroying pathogens, and also retards the gene WNT7A that is an immune regulatory gene involved in stem cell maintenance and healing when defective or mutated can slow the process of healing from infection. I extracted the EBVaNPCGOF genes manually as I stopped scraping sites long ago and this is a pub med article and would make no sense to build a program in a few days to scrape an article that could take a few hours at most with notetaking annotations and highlighting. The figure 5d, with the link analysis to the genes found in EBVaNPCGOF and in the GSEA study they found with Gene Ontology (GO) were of interest for my own research to use our database of differentially expressed genes we found and took the top and bottom genes with most upregulated and most down regulated changes from NPCto NPCwith EBV, as they used a commercial line of NPCand transfected the NPCwith their mutated LMP1 gene with H101R found in EBV. So we will do that in this part. I am leaving the previous work on the tail end of this part 2, usually I attach the subsequent parts the initial part of the project and separate with a bunch of equal signs to divide it or mark it with 3 or more astericks and I will do that on this project by dividing this part 2 from part 1 below with a bunch of equal signs and 3 astericks and label it part 1.
The genes found in the study of interest for this project to analyze and explore and possibly add or work with separately to see the changes or viability of them in the machine learning algorithm used are from the Figure 5D of the study. They also preferred using our same preferred supervised learning algorithm of random forest compared to others they used of gradient boosted machines, logistic regression, naive bayes, random forest, and a few unsupervised learning algorithms of support vector machines and K-Nearest Neighbors. The genes that we are interested in for this part 2 research analysis are the genes they found in EBVaNPCGOF and in their GSEA study with GO.
EBVaNPCGOF genes:
GSEA genes:
There were also some other genes mentioned in the study findings that weren’t included in the link analysis diagram of Figure 5d. Those genes are:
So, our list of genes to look up in our database before we extracted the most upregulated and most downregulated genes from NPCto NPCwith EBV will be: CD14, CSF1R, FOXP3, LTA, WNT7A, CD34, ANKRD11, IL7, LBH, PRRX1, HLA-DPA1, CEACAM6, and RAB38. Lets make this into a list of genes to extract.
studyGenes <- c("CD14", "CSF1R", "FOXP3", "LTA", "WNT7A", "CD34", "ANKRD11", "IL7", "LBH", "PRRX1", "HLA-DPA1", "CEACAM6", "RAB38")Lets read in that dataset and extract those genes into their own data set of top genes for that study GSE271486. We have to start from beginning as those data sets weren’t saved but no worries it was a quick analysis from start to end.
NPCEBV <- read.csv("GSE271486_HNE_1_MUT_LMP1_vs_HNE_1_WT_LMP1.csv", header=T, sep=',')
colnames(NPCEBV)## [1] "gene_id" "gene_name" "description"
## [4] "locus" "HNE_1_MUT_LMP1_1_count" "HNE_1_MUT_LMP1_2_count"
## [7] "HNE_1_MUT_LMP1_3_count" "HNE_1_WT_LMP1_1_count" "HNE_1_WT_LMP1_2_count"
## [10] "HNE_1_WT_LMP1_3_count" "HNE_1_MUT_LMP1_1_FPKM" "HNE_1_MUT_LMP1_2_FPKM"
## [13] "HNE_1_MUT_LMP1_3_FPKM" "HNE_1_WT_LMP1_1_FPKM" "HNE_1_WT_LMP1_2_FPKM"
## [16] "HNE_1_WT_LMP1_3_FPKM"NPCEBV$EBV_mean <- rowMeans(NPCEBV[,c(5:7)])
NPCEBV$baseline_mean <- rowMeans(NPCEBV[,c(8:10)])
NPCEBV$FoldchangeEBV_baseline <- NPCEBV$EBV_mean/NPCEBV$baseline_mean
NPCEBV_ordered <- NPCEBV[order(NPCEBV$FoldchangeEBV_baseline, decreasing=T),]
str(NPCEBV_ordered)## 'data.frame': 50868 obs. of 19 variables:
## $ gene_id : chr "ENSG00000001561" "ENSG00000004846" "ENSG00000005108" "ENSG00000006210" ...
## $ gene_name : chr "ENPP4" "ABCB5" "THSD7A" "CX3CL1" ...
## $ description : chr "ectonucleotide pyrophosphatase/phosphodiesterase 4 [Source:HGNC Symbol;Acc:HGNC:3359]" "ATP binding cassette subfamily B member 5 [Source:HGNC Symbol;Acc:HGNC:46]" "thrombospondin type 1 domain containing 7A [Source:HGNC Symbol;Acc:HGNC:22207]" "C-X3-C motif chemokine ligand 1 [Source:HGNC Symbol;Acc:HGNC:10647]" ...
## $ locus : chr "6:46129993-46146699:+" "7:20615207-20777038:+" "7:11370357-11832198:-" "16:57372458-57385048:+" ...
## $ HNE_1_MUT_LMP1_1_count: int 0 1 0 1 1 1 0 0 2 1 ...
## $ HNE_1_MUT_LMP1_2_count: int 1 0 0 0 0 0 0 2 0 0 ...
## $ HNE_1_MUT_LMP1_3_count: int 1 0 1 0 0 0 1 1 0 0 ...
## $ HNE_1_WT_LMP1_1_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_2_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_3_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_MUT_LMP1_1_FPKM : num 0 0.00529 0 0.00414 0.0065 ...
## $ HNE_1_MUT_LMP1_2_FPKM : num 0.0154 0 0 0 0 ...
## $ HNE_1_MUT_LMP1_3_FPKM : num 0.01476 0 0.00375 0 0 ...
## $ HNE_1_WT_LMP1_1_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_2_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_3_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ EBV_mean : num 0.667 0.333 0.333 0.333 0.333 ...
## $ baseline_mean : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FoldchangeEBV_baseline: num Inf Inf Inf Inf Inf ...Lets see if any of the study genes are in the gene_name feature of our data before we filtered it to exclude missing values or 0s in the samples and means.
common <- NPCEBV_ordered[which(NPCEBV_ordered$gene_name %in% studyGenes),]
common## gene_id gene_name
## 1699 ENSG00000086548 CEACAM6
## 669 ENSG00000049768 FOXP3
## 5431 ENSG00000123892 RAB38
## 12112 ENSG00000167522 ANKRD11
## 13644 ENSG00000174059 CD34
## 26608 ENSG00000231389 HLA-DPA1
## 12833 ENSG00000170458 CD14
## 15297 ENSG00000182578 CSF1R
## 19593 ENSG00000213626 LBH
## 2923 ENSG00000104432 IL7
## 4533 ENSG00000116132 PRRX1
## 9808 ENSG00000154764 WNT7A
## 23681 ENSG00000226979 LTA
## description
## 1699 CEA cell adhesion molecule 6 [Source:HGNC Symbol;Acc:HGNC:1818]
## 669 forkhead box P3 [Source:HGNC Symbol;Acc:HGNC:6106]
## 5431 RAB38, member RAS oncogene family [Source:HGNC Symbol;Acc:HGNC:9776]
## 12112 ankyrin repeat domain 11 [Source:HGNC Symbol;Acc:HGNC:21316]
## 13644 CD34 molecule [Source:HGNC Symbol;Acc:HGNC:1662]
## 26608 major histocompatibility complex, class II, DP alpha 1 [Source:HGNC Symbol;Acc:HGNC:4938]
## 12833 CD14 molecule [Source:HGNC Symbol;Acc:HGNC:1628]
## 15297 colony stimulating factor 1 receptor [Source:HGNC Symbol;Acc:HGNC:2433]
## 19593 LBH regulator of WNT signaling pathway [Source:HGNC Symbol;Acc:HGNC:29532]
## 2923 interleukin 7 [Source:HGNC Symbol;Acc:HGNC:6023]
## 4533 paired related homeobox 1 [Source:HGNC Symbol;Acc:HGNC:9142]
## 9808 Wnt family member 7A [Source:HGNC Symbol;Acc:HGNC:12786]
## 23681 lymphotoxin alpha [Source:HGNC Symbol;Acc:HGNC:6709]
## locus HNE_1_MUT_LMP1_1_count HNE_1_MUT_LMP1_2_count
## 1699 19:41750977-41772208:+ 3 3
## 669 X:49250436-49264826:- 6 7
## 5431 11:88113242-88175467:- 5 4
## 12112 16:89267627-89490561:- 6518 6902
## 13644 1:207880972-207911402:- 4 11
## 26608 6:33064569-33080775:- 0 1
## 12833 5:140631728-140633701:- 0 1
## 15297 5:150053291-150113372:- 3 2
## 19593 2:30231531-30323730:+ 8 10
## 2923 8:78675743-78805523:- 0 0
## 4533 1:170662728-170739419:+ 0 0
## 9808 3:13816258-13880121:- 0 0
## 23681 6:31572054-31574324:+ 0 0
## HNE_1_MUT_LMP1_3_count HNE_1_WT_LMP1_1_count HNE_1_WT_LMP1_2_count
## 1699 2 0 0
## 669 4 3 1
## 5431 3 3 2
## 12112 3838 3197 2469
## 13644 5 13 10
## 26608 2 2 3
## 12833 3 4 3
## 15297 4 16 7
## 19593 10 50 19
## 2923 1 2 1
## 4533 1 2 2
## 9808 0 2 1
## 23681 0 2 1
## HNE_1_WT_LMP1_3_count HNE_1_MUT_LMP1_1_FPKM HNE_1_MUT_LMP1_2_FPKM
## 1699 0 0.04930920 0.065840603
## 669 2 0.06943194 0.108171793
## 5431 0 0.08618129 0.092081288
## 12112 3724 15.22033021 21.523280160
## 13644 17 0.01727330 0.063439477
## 26608 3 0.00000000 0.007629582
## 12833 5 0.00000000 0.018610955
## 15297 6 0.02419598 0.021542269
## 19593 26 0.05776895 0.096431264
## 2923 2 0.00000000 0.000000000
## 4533 4 0.00000000 0.000000000
## 9808 2 0.00000000 0.000000000
## 23681 3 0.00000000 0.000000000
## HNE_1_MUT_LMP1_3_FPKM HNE_1_WT_LMP1_1_FPKM HNE_1_WT_LMP1_2_FPKM
## 1699 0.041948267 0.000000000 0.000000000
## 669 0.059082292 0.025135938 0.011115633
## 5431 0.066000178 0.037441990 0.033115249
## 12112 11.439281190 5.405314676 5.537387379
## 13644 0.027560037 0.040650254 0.041479696
## 26608 0.014587793 0.008277097 0.016466422
## 12833 0.053385517 0.040380901 0.040166796
## 15297 0.041178567 0.093428733 0.054216856
## 19593 0.092169972 0.261410758 0.131771024
## 2923 0.009805968 0.011122054 0.007377604
## 4533 0.006280126 0.007122999 0.009449813
## 9808 0.000000000 0.015655415 0.010384724
## 23681 0.000000000 0.036093823 0.023942156
## HNE_1_WT_LMP1_3_FPKM EBV_mean baseline_mean FoldchangeEBV_baseline
## 1699 0.00000000 2.6666667 0.000000 Inf
## 669 0.02852506 5.6666667 2.000000 2.8333333
## 5431 0.00000000 4.0000000 1.666667 2.4000000
## 12112 10.71962836 5752.6666667 3130.000000 1.8379127
## 13644 0.09050085 6.6666667 13.333333 0.5000000
## 26608 0.02113450 1.0000000 2.666667 0.3750000
## 12833 0.08592286 1.3333333 4.000000 0.3333333
## 15297 0.05964330 3.0000000 9.666667 0.3103448
## 19593 0.23143041 9.3333333 31.666667 0.2947368
## 2923 0.01893249 0.3333333 1.666667 0.2000000
## 4533 0.02425021 0.3333333 2.666667 0.1250000
## 9808 0.02664939 0.0000000 1.666667 0.0000000
## 23681 0.09216093 0.0000000 2.000000 0.0000000All 13 genes are in the data we have before filtering it. Lets see if the 13 genes are also in the data once we excluded the 0s for mean values and fold change irregularities created by having 0s for the baseline or EBV mutated mean value of genes.
NPCEBV2 <- subset(NPCEBV_ordered,NPCEBV_ordered$baseline_mean > 0 & NPCEBV_ordered$EBV_mean > 0)
str(NPCEBV2)## 'data.frame': 21932 obs. of 19 variables:
## $ gene_id : chr "ENSG00000262660" "ENSG00000267303" "ENSG00000198211" "ENSG00000240963" ...
## $ gene_name : chr "AC139530.2" "AC011511.4" "AC092143.1" "AL645465.1" ...
## $ description : chr "novel protein" "novel transcript" "novel protein (MC1R-TUBB3 readthrough)" "novel transcript, antisense to C1orf100" ...
## $ locus : chr "17:81703371-81720539:+" "19:10315471-10320678:-" "16:89919165-89936092:+" "1:244375100-244409592:-" ...
## $ HNE_1_MUT_LMP1_1_count: int 520 617 211 93 1 15 0 0 5 0 ...
## $ HNE_1_MUT_LMP1_2_count: int 0 0 0 3 63 203 36 54 4 13 ...
## $ HNE_1_MUT_LMP1_3_count: int 0 0 62 85 49 12 32 0 7 1 ...
## $ HNE_1_WT_LMP1_1_count : int 0 0 3 1 2 4 0 3 0 0 ...
## $ HNE_1_WT_LMP1_2_count : int 1 0 0 1 1 2 3 0 0 0 ...
## $ HNE_1_WT_LMP1_3_count : int 0 4 0 2 1 3 0 0 1 1 ...
## $ HNE_1_MUT_LMP1_1_FPKM : num 14.5446 17.9294 4.0862 7.7876 0.0913 ...
## $ HNE_1_MUT_LMP1_2_FPKM : num 0 0 0 0.335 7.679 ...
## $ HNE_1_MUT_LMP1_3_FPKM : num 0 0 1.53 9.08 5.71 ...
## $ HNE_1_WT_LMP1_1_FPKM : num 0 0 0.0421 0.0606 0.1322 ...
## $ HNE_1_WT_LMP1_2_FPKM : num 0.0269 0 0 0.0804 0.0877 ...
## $ HNE_1_WT_LMP1_3_FPKM : num 0 0.143 0 0.206 0.113 ...
## $ EBV_mean : num 173.3 205.7 91 60.3 37.7 ...
## $ baseline_mean : num 0.333 1.333 1 1.333 1.333 ...
## $ FoldchangeEBV_baseline: num 520 154.2 91 45.3 28.2 ...range(NPCEBV2$FoldchangeEBV_baseline)## [1] 0.01107011 520.00000000commonFiltered <- NPCEBV2[which(NPCEBV2$gene_name %in% studyGenes),]
commonFiltered## gene_id gene_name
## 669 ENSG00000049768 FOXP3
## 5431 ENSG00000123892 RAB38
## 12112 ENSG00000167522 ANKRD11
## 13644 ENSG00000174059 CD34
## 26608 ENSG00000231389 HLA-DPA1
## 12833 ENSG00000170458 CD14
## 15297 ENSG00000182578 CSF1R
## 19593 ENSG00000213626 LBH
## 2923 ENSG00000104432 IL7
## 4533 ENSG00000116132 PRRX1
## description
## 669 forkhead box P3 [Source:HGNC Symbol;Acc:HGNC:6106]
## 5431 RAB38, member RAS oncogene family [Source:HGNC Symbol;Acc:HGNC:9776]
## 12112 ankyrin repeat domain 11 [Source:HGNC Symbol;Acc:HGNC:21316]
## 13644 CD34 molecule [Source:HGNC Symbol;Acc:HGNC:1662]
## 26608 major histocompatibility complex, class II, DP alpha 1 [Source:HGNC Symbol;Acc:HGNC:4938]
## 12833 CD14 molecule [Source:HGNC Symbol;Acc:HGNC:1628]
## 15297 colony stimulating factor 1 receptor [Source:HGNC Symbol;Acc:HGNC:2433]
## 19593 LBH regulator of WNT signaling pathway [Source:HGNC Symbol;Acc:HGNC:29532]
## 2923 interleukin 7 [Source:HGNC Symbol;Acc:HGNC:6023]
## 4533 paired related homeobox 1 [Source:HGNC Symbol;Acc:HGNC:9142]
## locus HNE_1_MUT_LMP1_1_count HNE_1_MUT_LMP1_2_count
## 669 X:49250436-49264826:- 6 7
## 5431 11:88113242-88175467:- 5 4
## 12112 16:89267627-89490561:- 6518 6902
## 13644 1:207880972-207911402:- 4 11
## 26608 6:33064569-33080775:- 0 1
## 12833 5:140631728-140633701:- 0 1
## 15297 5:150053291-150113372:- 3 2
## 19593 2:30231531-30323730:+ 8 10
## 2923 8:78675743-78805523:- 0 0
## 4533 1:170662728-170739419:+ 0 0
## HNE_1_MUT_LMP1_3_count HNE_1_WT_LMP1_1_count HNE_1_WT_LMP1_2_count
## 669 4 3 1
## 5431 3 3 2
## 12112 3838 3197 2469
## 13644 5 13 10
## 26608 2 2 3
## 12833 3 4 3
## 15297 4 16 7
## 19593 10 50 19
## 2923 1 2 1
## 4533 1 2 2
## HNE_1_WT_LMP1_3_count HNE_1_MUT_LMP1_1_FPKM HNE_1_MUT_LMP1_2_FPKM
## 669 2 0.06943194 0.108171793
## 5431 0 0.08618129 0.092081288
## 12112 3724 15.22033021 21.523280160
## 13644 17 0.01727330 0.063439477
## 26608 3 0.00000000 0.007629582
## 12833 5 0.00000000 0.018610955
## 15297 6 0.02419598 0.021542269
## 19593 26 0.05776895 0.096431264
## 2923 2 0.00000000 0.000000000
## 4533 4 0.00000000 0.000000000
## HNE_1_MUT_LMP1_3_FPKM HNE_1_WT_LMP1_1_FPKM HNE_1_WT_LMP1_2_FPKM
## 669 0.059082292 0.025135938 0.011115633
## 5431 0.066000178 0.037441990 0.033115249
## 12112 11.439281190 5.405314676 5.537387379
## 13644 0.027560037 0.040650254 0.041479696
## 26608 0.014587793 0.008277097 0.016466422
## 12833 0.053385517 0.040380901 0.040166796
## 15297 0.041178567 0.093428733 0.054216856
## 19593 0.092169972 0.261410758 0.131771024
## 2923 0.009805968 0.011122054 0.007377604
## 4533 0.006280126 0.007122999 0.009449813
## HNE_1_WT_LMP1_3_FPKM EBV_mean baseline_mean FoldchangeEBV_baseline
## 669 0.02852506 5.6666667 2.000000 2.8333333
## 5431 0.00000000 4.0000000 1.666667 2.4000000
## 12112 10.71962836 5752.6666667 3130.000000 1.8379127
## 13644 0.09050085 6.6666667 13.333333 0.5000000
## 26608 0.02113450 1.0000000 2.666667 0.3750000
## 12833 0.08592286 1.3333333 4.000000 0.3333333
## 15297 0.05964330 3.0000000 9.666667 0.3103448
## 19593 0.23143041 9.3333333 31.666667 0.2947368
## 2923 0.01893249 0.3333333 1.666667 0.2000000
## 4533 0.02425021 0.3333333 2.666667 0.1250000Only 10 of the genes are in the data once 0s were removed for genes that created irregularities in the fold change values of EBV mutated NCP compared to NCP only.
Lets see those genes and which removed to compare information.
c <- common$gene_name
c## [1] "CEACAM6" "FOXP3" "RAB38" "ANKRD11" "CD34" "HLA-DPA1"
## [7] "CD14" "CSF1R" "LBH" "IL7" "PRRX1" "WNT7A"
## [13] "LTA"cf <- commonFiltered$gene_name
cf## [1] "FOXP3" "RAB38" "ANKRD11" "CD34" "HLA-DPA1" "CD14"
## [7] "CSF1R" "LBH" "IL7" "PRRX1"nc <- c[-which(c %in% cf)]
nc## [1] "CEACAM6" "WNT7A" "LTA"There are 3 genes that aren’t in the study gene set once filtered to exclude 0 fold change values or irregular fold change values from having 0 mean value for the gene. Lets look at the data of these 3 genes from unfiltered data.
genes3 <- NPCEBV_ordered[which(NPCEBV_ordered$gene_name %in% nc),]
genes3## gene_id gene_name
## 1699 ENSG00000086548 CEACAM6
## 9808 ENSG00000154764 WNT7A
## 23681 ENSG00000226979 LTA
## description
## 1699 CEA cell adhesion molecule 6 [Source:HGNC Symbol;Acc:HGNC:1818]
## 9808 Wnt family member 7A [Source:HGNC Symbol;Acc:HGNC:12786]
## 23681 lymphotoxin alpha [Source:HGNC Symbol;Acc:HGNC:6709]
## locus HNE_1_MUT_LMP1_1_count HNE_1_MUT_LMP1_2_count
## 1699 19:41750977-41772208:+ 3 3
## 9808 3:13816258-13880121:- 0 0
## 23681 6:31572054-31574324:+ 0 0
## HNE_1_MUT_LMP1_3_count HNE_1_WT_LMP1_1_count HNE_1_WT_LMP1_2_count
## 1699 2 0 0
## 9808 0 2 1
## 23681 0 2 1
## HNE_1_WT_LMP1_3_count HNE_1_MUT_LMP1_1_FPKM HNE_1_MUT_LMP1_2_FPKM
## 1699 0 0.0493092 0.0658406
## 9808 2 0.0000000 0.0000000
## 23681 3 0.0000000 0.0000000
## HNE_1_MUT_LMP1_3_FPKM HNE_1_WT_LMP1_1_FPKM HNE_1_WT_LMP1_2_FPKM
## 1699 0.04194827 0.00000000 0.00000000
## 9808 0.00000000 0.01565541 0.01038472
## 23681 0.00000000 0.03609382 0.02394216
## HNE_1_WT_LMP1_3_FPKM EBV_mean baseline_mean FoldchangeEBV_baseline
## 1699 0.00000000 2.666667 0.000000 Inf
## 9808 0.02664939 0.000000 1.666667 0
## 23681 0.09216093 0.000000 2.000000 0Looks like the 2 key genes of CEACAM6 and WNT7A didn’t make the filter as well as LTA, because CEACAM6 was highly upregulated from almost 0 (we can say almost 0 because in gene expression data after normalizing the data between 0 and 1 some values can be 1X1.0^-12 or 0.000000000001) and excluded or rounded to 0.00 automatically by the computer that rounds to the 8th or 12th figure). The WNT7A gene is found to be highly downregulated in EBV mutated NCP as the study claimed and it is responsible for stem cell maintenance (WBCs, RBCs, and lymphocytes are stem cells and natural killer cells are WBCs of monocytes, eosinophiles, and neutrophiles) and renewal, while CEACAM6 the study cited a source as saying it is upregulated in many cancers and proliferation and metastasis of tumors we can see here is highly upregulated. And LT7 we can go to the included summary of the gene in the gene data they provided that it is highly downregulated and is a lymphotoxin alpha which sounds like that is a good thing for overall immunity. Lets look at all the genes next to their fold change with only the gene name, the description, and the fold change value.
studyGenesDF <- NPCEBV_ordered[which(NPCEBV_ordered$gene_name %in% studyGenes), c(2,3,19)]
studyGenesDF## gene_name
## 1699 CEACAM6
## 669 FOXP3
## 5431 RAB38
## 12112 ANKRD11
## 13644 CD34
## 26608 HLA-DPA1
## 12833 CD14
## 15297 CSF1R
## 19593 LBH
## 2923 IL7
## 4533 PRRX1
## 9808 WNT7A
## 23681 LTA
## description
## 1699 CEA cell adhesion molecule 6 [Source:HGNC Symbol;Acc:HGNC:1818]
## 669 forkhead box P3 [Source:HGNC Symbol;Acc:HGNC:6106]
## 5431 RAB38, member RAS oncogene family [Source:HGNC Symbol;Acc:HGNC:9776]
## 12112 ankyrin repeat domain 11 [Source:HGNC Symbol;Acc:HGNC:21316]
## 13644 CD34 molecule [Source:HGNC Symbol;Acc:HGNC:1662]
## 26608 major histocompatibility complex, class II, DP alpha 1 [Source:HGNC Symbol;Acc:HGNC:4938]
## 12833 CD14 molecule [Source:HGNC Symbol;Acc:HGNC:1628]
## 15297 colony stimulating factor 1 receptor [Source:HGNC Symbol;Acc:HGNC:2433]
## 19593 LBH regulator of WNT signaling pathway [Source:HGNC Symbol;Acc:HGNC:29532]
## 2923 interleukin 7 [Source:HGNC Symbol;Acc:HGNC:6023]
## 4533 paired related homeobox 1 [Source:HGNC Symbol;Acc:HGNC:9142]
## 9808 Wnt family member 7A [Source:HGNC Symbol;Acc:HGNC:12786]
## 23681 lymphotoxin alpha [Source:HGNC Symbol;Acc:HGNC:6709]
## FoldchangeEBV_baseline
## 1699 Inf
## 669 2.8333333
## 5431 2.4000000
## 12112 1.8379127
## 13644 0.5000000
## 26608 0.3750000
## 12833 0.3333333
## 15297 0.3103448
## 19593 0.2947368
## 2923 0.2000000
## 4533 0.1250000
## 9808 0.0000000
## 23681 0.0000000The top 4 genes are upregulated while the last 9 genes are downregulated in comparing the EBV mutated H101R of LMP1 gene in NPC to NPC without the mutation or the wild type. The top 4 upregulated genes includes CEACAM6, FOXP3, and RAB38 as well as ANKRD11. CEACAM6 is seen in tumor creation, proliferation, and metastasis, while FOXP3 is seen in T cell anergy or lazy non proliferative T cells like CD4 that gets Natural Killer Cells (NKC) to fight infection and destroy contaminants, and that RAB38 that is seen highly upregulated in the phenotypic metastatic pancreatic cancer is upregulated but only at a value of 2.4 which could be almost 2 1/2 times the gene expression it should be at in controlled NPC.
Lets write both of these data sets to the csv format and upload and share to view or use at a later time.
write.csv(NPCEBV_ordered, "NPCEBV_ordered.csv", row.names=F)
write.csv(studyGenesDF, "studyGenesDF.csv",row.names=F)
write.csv(NPCEBV2, "NCPEBV_filtered.csv", row.names=F)You can get the ordered and not filtered data set of 50,868 genes here, the data set of only the 13 genes in the study from unfiltered data set here, and the filtered data of genes that is 21,932 genes long here.
Thanks so much, beyond this we could test to see how well these study genes perform in predicting the EBVaNCP or the NCP class samples.
================================================
*** Part 1
In this project, we seek to understand target genes in the study GSE271486 that is a recently published study analyzing nasopharyngeal carcinoma (NCP) and the effects of Epstein-Barr Viral (EBV) infection to changes in the latent membrane protein 1 or LMP1. The study used machine learning to predict a model developed from 26 EBV infected patients that showed there were changes in oncogenes that could help indicate NPC developing in those infected with EBV. Their gene targets found through the gene ontology fingerprint of EBV and NPC were associated with changes in H101R and its affects on immune regulation processes to the FOXP3-T cell anergy and WNT7A-stem cell population maintenance.
This study, GSE271486, essentially found the most differentially expressed genes in an infected EBV population of 26, then used those gene targets and gene groups from the Gene Ontology fingerprint of NPC and EBV infection to determine the top feature in identifying EBV infection was H101R mutation or genes of most variability as in cell differentiation in cancer indicator using a commercial cell line of NCP called HNE-1 that they infected with the mutant H101R in vitro or in a petri dish and watched changes. They compared the non-infected NPC to the EBV infected NCP petri dish outcomes after 7 days, and confirmed their findings that mutations in H101R to the LMBP-1 region gave a higher risk of NCP. However, the cell line alread had NCP, so we can look at the data they supplied for us already in a csv format that has 3 samples of infected EBV called HNE-1MUT-LMP1 and 3 not infected assumed as it is the other group HNE-1WT_LMP1. There are 3 counts and 3 fragments per kilomillion features (fpkm) for each group of infected EBV or not infected. I expect to see the group of most differentially expressed genes being those related to EBV and NCP fingerprint gene groups from the Gene Ontology Fingerprint.
Sounds like a super fascinating study, lets see what we find and if we can use these top genes in our database of pathologies that we can use to build a machine model to predict EBV associated pathologies, EBV infection, fibromyalgia, Lyme disease, multiple sclerosis, mononucleosis, Hodgkin’s lymphoma, Burkett’s lymphoma, Large B-cell lymphoma, or Natural Killer T-cell Lymphoma.
Lets read in the csv file that was provided.
ncpEBV <- read.csv("GSE271486_HNE_1_MUT_LMP1_vs_HNE_1_WT_LMP1.csv", header=T, sep=',')
colnames(ncpEBV)## [1] "gene_id" "gene_name" "description"
## [4] "locus" "HNE_1_MUT_LMP1_1_count" "HNE_1_MUT_LMP1_2_count"
## [7] "HNE_1_MUT_LMP1_3_count" "HNE_1_WT_LMP1_1_count" "HNE_1_WT_LMP1_2_count"
## [10] "HNE_1_WT_LMP1_3_count" "HNE_1_MUT_LMP1_1_FPKM" "HNE_1_MUT_LMP1_2_FPKM"
## [13] "HNE_1_MUT_LMP1_3_FPKM" "HNE_1_WT_LMP1_1_FPKM" "HNE_1_WT_LMP1_2_FPKM"
## [16] "HNE_1_WT_LMP1_3_FPKM"head(ncpEBV)## gene_id gene_name
## 1 ENSG00000000003 TSPAN6
## 2 ENSG00000000005 TNMD
## 3 ENSG00000000419 DPM1
## 4 ENSG00000000457 SCYL3
## 5 ENSG00000000460 C1orf112
## 6 ENSG00000000938 FGR
## description
## 1 tetraspanin 6 [Source:HGNC Symbol;Acc:HGNC:11858]
## 2 tenomodulin [Source:HGNC Symbol;Acc:HGNC:17757]
## 3 dolichyl-phosphate mannosyltransferase subunit 1, catalytic [Source:HGNC Symbol;Acc:HGNC:3005]
## 4 SCY1 like pseudokinase 3 [Source:HGNC Symbol;Acc:HGNC:19285]
## 5 chromosome 1 open reading frame 112 [Source:HGNC Symbol;Acc:HGNC:25565]
## 6 FGR proto-oncogene, Src family tyrosine kinase [Source:HGNC Symbol;Acc:HGNC:3697]
## locus HNE_1_MUT_LMP1_1_count HNE_1_MUT_LMP1_2_count
## 1 X:100627109-100639991:- 1057 700
## 2 X:100584802-100599885:+ 0 0
## 3 20:50934867-50958555:- 937 678
## 4 1:169849631-169894267:- 200 136
## 5 1:169662007-169854080:+ 537 474
## 6 1:27612064-27635277:- 0 0
## HNE_1_MUT_LMP1_3_count HNE_1_WT_LMP1_1_count HNE_1_WT_LMP1_2_count
## 1 855 1643 826
## 2 0 0 0
## 3 932 1114 930
## 4 158 273 120
## 5 471 839 617
## 6 0 0 0
## HNE_1_WT_LMP1_3_count HNE_1_MUT_LMP1_1_FPKM HNE_1_MUT_LMP1_2_FPKM
## 1 875 9.355201 8.273695
## 2 0 0.000000 0.000000
## 3 845 19.973180 19.300172
## 4 177 1.247310 1.132686
## 5 594 3.772203 4.446549
## 6 0 0.000000 0.000000
## HNE_1_MUT_LMP1_3_FPKM HNE_1_WT_LMP1_1_FPKM HNE_1_WT_LMP1_2_FPKM
## 1 9.658906 10.528930 7.0215334
## 2 0.000000 0.000000 0.0000000
## 3 25.357562 17.193387 19.0398968
## 4 1.257724 1.232754 0.7187855
## 5 4.223044 4.267284 4.1627526
## 6 0.000000 0.000000 0.0000000
## HNE_1_WT_LMP1_3_FPKM
## 1 9.546551
## 2 0.000000
## 3 22.203657
## 4 1.360751
## 5 5.143604
## 6 0.000000str(ncpEBV)## 'data.frame': 50868 obs. of 16 variables:
## $ gene_id : chr "ENSG00000000003" "ENSG00000000005" "ENSG00000000419" "ENSG00000000457" ...
## $ gene_name : chr "TSPAN6" "TNMD" "DPM1" "SCYL3" ...
## $ description : chr "tetraspanin 6 [Source:HGNC Symbol;Acc:HGNC:11858]" "tenomodulin [Source:HGNC Symbol;Acc:HGNC:17757]" "dolichyl-phosphate mannosyltransferase subunit 1, catalytic [Source:HGNC Symbol;Acc:HGNC:3005]" "SCY1 like pseudokinase 3 [Source:HGNC Symbol;Acc:HGNC:19285]" ...
## $ locus : chr "X:100627109-100639991:-" "X:100584802-100599885:+" "20:50934867-50958555:-" "1:169849631-169894267:-" ...
## $ HNE_1_MUT_LMP1_1_count: int 1057 0 937 200 537 0 277 2475 1635 796 ...
## $ HNE_1_MUT_LMP1_2_count: int 700 0 678 136 474 0 256 1634 1753 668 ...
## $ HNE_1_MUT_LMP1_3_count: int 855 0 932 158 471 0 316 1864 1833 804 ...
## $ HNE_1_WT_LMP1_1_count : int 1643 0 1114 273 839 0 386 3631 2044 1101 ...
## $ HNE_1_WT_LMP1_2_count : int 826 0 930 120 617 0 299 1907 1732 878 ...
## $ HNE_1_WT_LMP1_3_count : int 875 0 845 177 594 0 280 2069 2159 774 ...
## $ HNE_1_MUT_LMP1_1_FPKM : num 9.36 0 19.97 1.25 3.77 ...
## $ HNE_1_MUT_LMP1_2_FPKM : num 8.27 0 19.3 1.13 4.45 ...
## $ HNE_1_MUT_LMP1_3_FPKM : num 9.66 0 25.36 1.26 4.22 ...
## $ HNE_1_WT_LMP1_1_FPKM : num 10.53 0 17.19 1.23 4.27 ...
## $ HNE_1_WT_LMP1_2_FPKM : num 7.022 0 19.04 0.719 4.163 ...
## $ HNE_1_WT_LMP1_3_FPKM : num 9.55 0 22.2 1.36 5.14 ...Lets do a quick analysis with fold change values and see the most differentially expressed genes using only the counts (no filtering with counts and FPKM or mitochondrial DNA percent as was done in data made for Seurat) from non infected EBV to infected EBV on the NCP commercial cell line.
ncpEBV$EBV_mean <- rowMeans(ncpEBV[,c(5:7)])
ncpEBV$baseline_mean <- rowMeans(ncpEBV[,c(8:10)])
ncpEBV$FoldchangeEBV_baseline <- ncpEBV$EBV_mean/ncpEBV$baseline_mean
ncpEBV_ordered <- ncpEBV[order(ncpEBV$FoldchangeEBV_baseline, decreasing=T),]
str(ncpEBV_ordered)## 'data.frame': 50868 obs. of 19 variables:
## $ gene_id : chr "ENSG00000001561" "ENSG00000004846" "ENSG00000005108" "ENSG00000006210" ...
## $ gene_name : chr "ENPP4" "ABCB5" "THSD7A" "CX3CL1" ...
## $ description : chr "ectonucleotide pyrophosphatase/phosphodiesterase 4 [Source:HGNC Symbol;Acc:HGNC:3359]" "ATP binding cassette subfamily B member 5 [Source:HGNC Symbol;Acc:HGNC:46]" "thrombospondin type 1 domain containing 7A [Source:HGNC Symbol;Acc:HGNC:22207]" "C-X3-C motif chemokine ligand 1 [Source:HGNC Symbol;Acc:HGNC:10647]" ...
## $ locus : chr "6:46129993-46146699:+" "7:20615207-20777038:+" "7:11370357-11832198:-" "16:57372458-57385048:+" ...
## $ HNE_1_MUT_LMP1_1_count: int 0 1 0 1 1 1 0 0 2 1 ...
## $ HNE_1_MUT_LMP1_2_count: int 1 0 0 0 0 0 0 2 0 0 ...
## $ HNE_1_MUT_LMP1_3_count: int 1 0 1 0 0 0 1 1 0 0 ...
## $ HNE_1_WT_LMP1_1_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_2_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_3_count : int 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_MUT_LMP1_1_FPKM : num 0 0.00529 0 0.00414 0.0065 ...
## $ HNE_1_MUT_LMP1_2_FPKM : num 0.0154 0 0 0 0 ...
## $ HNE_1_MUT_LMP1_3_FPKM : num 0.01476 0 0.00375 0 0 ...
## $ HNE_1_WT_LMP1_1_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_2_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ HNE_1_WT_LMP1_3_FPKM : num 0 0 0 0 0 0 0 0 0 0 ...
## $ EBV_mean : num 0.667 0.333 0.333 0.333 0.333 ...
## $ baseline_mean : num 0 0 0 0 0 0 0 0 0 0 ...
## $ FoldchangeEBV_baseline: num Inf Inf Inf Inf Inf ...ncpEBV2 <- subset(ncpEBV_ordered,ncpEBV_ordered$baseline_mean > 0 & ncpEBV_ordered$EBV_mean > 0)
str(ncpEBV2)## 'data.frame': 21932 obs. of 19 variables:
## $ gene_id : chr "ENSG00000262660" "ENSG00000267303" "ENSG00000198211" "ENSG00000240963" ...
## $ gene_name : chr "AC139530.2" "AC011511.4" "AC092143.1" "AL645465.1" ...
## $ description : chr "novel protein" "novel transcript" "novel protein (MC1R-TUBB3 readthrough)" "novel transcript, antisense to C1orf100" ...
## $ locus : chr "17:81703371-81720539:+" "19:10315471-10320678:-" "16:89919165-89936092:+" "1:244375100-244409592:-" ...
## $ HNE_1_MUT_LMP1_1_count: int 520 617 211 93 1 15 0 0 5 0 ...
## $ HNE_1_MUT_LMP1_2_count: int 0 0 0 3 63 203 36 54 4 13 ...
## $ HNE_1_MUT_LMP1_3_count: int 0 0 62 85 49 12 32 0 7 1 ...
## $ HNE_1_WT_LMP1_1_count : int 0 0 3 1 2 4 0 3 0 0 ...
## $ HNE_1_WT_LMP1_2_count : int 1 0 0 1 1 2 3 0 0 0 ...
## $ HNE_1_WT_LMP1_3_count : int 0 4 0 2 1 3 0 0 1 1 ...
## $ HNE_1_MUT_LMP1_1_FPKM : num 14.5446 17.9294 4.0862 7.7876 0.0913 ...
## $ HNE_1_MUT_LMP1_2_FPKM : num 0 0 0 0.335 7.679 ...
## $ HNE_1_MUT_LMP1_3_FPKM : num 0 0 1.53 9.08 5.71 ...
## $ HNE_1_WT_LMP1_1_FPKM : num 0 0 0.0421 0.0606 0.1322 ...
## $ HNE_1_WT_LMP1_2_FPKM : num 0.0269 0 0 0.0804 0.0877 ...
## $ HNE_1_WT_LMP1_3_FPKM : num 0 0.143 0 0.206 0.113 ...
## $ EBV_mean : num 173.3 205.7 91 60.3 37.7 ...
## $ baseline_mean : num 0.333 1.333 1 1.333 1.333 ...
## $ FoldchangeEBV_baseline: num 520 154.2 91 45.3 28.2 ...We have a dataset that is smaller than original by excluding those entries where the mean of EBV infected and the mean of the baseline NCP samples are higher than 0.
range(ncpEBV2$FoldchangeEBV_baseline)## [1] 0.01107011 520.00000000The range is all positive but those less than 1 are underexpressed genes or inhibited in EBV infected NCP compared to baseline NCP samples. Lets see what are top 10 overexpressed genes are in stimulated the most or inhibited the most for top and bottom 10 genes. We should also look at a list of those exact genes that the study used as the ‘Gene Ontology Fingerprint’ for ‘EBVaNPC GOF’ perhaps with an internet search. The internet returned the Pubmed article by this GSE study.
topBottom10stimulatedInhibited <- ncpEBV2[c(1:10,21923:21932),]
topBottom10stimulatedInhibited[,c(2,19)]## gene_name FoldchangeEBV_baseline
## 42404 AC139530.2 520.00000000
## 43879 AC011511.4 154.25000000
## 17708 AC092143.1 91.00000000
## 31972 AL645465.1 45.25000000
## 43987 AL353997.3 28.25000000
## 46317 AC022384.1 25.55555556
## 95 MPO 22.66666667
## 39685 TRIM6-TRIM34 18.00000000
## 17107 FUT4 16.00000000
## 4383 ITGA4 14.00000000
## 23008 RPL23AP34 0.07692308
## 40275 AC027088.3 0.07692308
## 8269 MYOM3 0.07142857
## 70 CALCR 0.06666667
## 1480 MOXD1 0.06250000
## 11402 SMAD5-AS1 0.05882353
## 45524 HSPE1P7 0.04761905
## 40946 AC026464.2 0.03557312
## 12885 RASA4B 0.03225806
## 42231 AC026954.2 0.01107011Lets make a matrix and split into testing and training sets comparing the samples, this will be small since only 3 baseline and 3 EBV infected NCP.
colnames(topBottom10stimulatedInhibited)## [1] "gene_id" "gene_name" "description"
## [4] "locus" "HNE_1_MUT_LMP1_1_count" "HNE_1_MUT_LMP1_2_count"
## [7] "HNE_1_MUT_LMP1_3_count" "HNE_1_WT_LMP1_1_count" "HNE_1_WT_LMP1_2_count"
## [10] "HNE_1_WT_LMP1_3_count" "HNE_1_MUT_LMP1_1_FPKM" "HNE_1_MUT_LMP1_2_FPKM"
## [13] "HNE_1_MUT_LMP1_3_FPKM" "HNE_1_WT_LMP1_1_FPKM" "HNE_1_WT_LMP1_2_FPKM"
## [16] "HNE_1_WT_LMP1_3_FPKM" "EBV_mean" "baseline_mean"
## [19] "FoldchangeEBV_baseline"header <- topBottom10stimulatedInhibited$gene_name
class <- c('EBVncp','EBVncp','EBVncp','ncp','ncp','ncp')
df <- topBottom10stimulatedInhibited[,c(5:10)]
str(df)## 'data.frame': 20 obs. of 6 variables:
## $ HNE_1_MUT_LMP1_1_count: int 520 617 211 93 1 15 0 0 5 0 ...
## $ HNE_1_MUT_LMP1_2_count: int 0 0 0 3 63 203 36 54 4 13 ...
## $ HNE_1_MUT_LMP1_3_count: int 0 0 62 85 49 12 32 0 7 1 ...
## $ HNE_1_WT_LMP1_1_count : int 0 0 3 1 2 4 0 3 0 0 ...
## $ HNE_1_WT_LMP1_2_count : int 1 0 0 1 1 2 3 0 0 0 ...
## $ HNE_1_WT_LMP1_3_count : int 0 4 0 2 1 3 0 0 1 1 ...ncpEBV_matrix <- data.frame(t(df))
colnames(ncpEBV_matrix) <- header
ncpEBV_matrix$class <- class
str(ncpEBV_matrix)## 'data.frame': 6 obs. of 21 variables:
## $ AC139530.2 : int 520 0 0 0 1 0
## $ AC011511.4 : int 617 0 0 0 0 4
## $ AC092143.1 : int 211 0 62 3 0 0
## $ AL645465.1 : int 93 3 85 1 1 2
## $ AL353997.3 : int 1 63 49 2 1 1
## $ AC022384.1 : int 15 203 12 4 2 3
## $ MPO : int 0 36 32 0 3 0
## $ TRIM6-TRIM34: int 0 54 0 3 0 0
## $ FUT4 : int 5 4 7 0 0 1
## $ ITGA4 : int 0 13 1 0 0 1
## $ RPL23AP34 : int 0 1 0 5 1 7
## $ AC027088.3 : int 2 0 0 9 13 4
## $ MYOM3 : int 1 0 0 7 7 0
## $ CALCR : int 0 0 1 6 1 8
## $ MOXD1 : int 1 0 0 12 2 2
## $ SMAD5-AS1 : int 0 1 0 14 2 1
## $ HSPE1P7 : int 0 1 0 0 15 6
## $ AC026464.2 : int 4 2 3 184 3 66
## $ RASA4B : int 2 0 0 59 3 0
## $ AC026954.2 : int 1 1 1 269 1 1
## $ class : chr "EBVncp" "EBVncp" "EBVncp" "ncp" ...ncpEBV_matrix$class <- as.factor(ncpEBV_matrix$class)
str(ncpEBV_matrix)## 'data.frame': 6 obs. of 21 variables:
## $ AC139530.2 : int 520 0 0 0 1 0
## $ AC011511.4 : int 617 0 0 0 0 4
## $ AC092143.1 : int 211 0 62 3 0 0
## $ AL645465.1 : int 93 3 85 1 1 2
## $ AL353997.3 : int 1 63 49 2 1 1
## $ AC022384.1 : int 15 203 12 4 2 3
## $ MPO : int 0 36 32 0 3 0
## $ TRIM6-TRIM34: int 0 54 0 3 0 0
## $ FUT4 : int 5 4 7 0 0 1
## $ ITGA4 : int 0 13 1 0 0 1
## $ RPL23AP34 : int 0 1 0 5 1 7
## $ AC027088.3 : int 2 0 0 9 13 4
## $ MYOM3 : int 1 0 0 7 7 0
## $ CALCR : int 0 0 1 6 1 8
## $ MOXD1 : int 1 0 0 12 2 2
## $ SMAD5-AS1 : int 0 1 0 14 2 1
## $ HSPE1P7 : int 0 1 0 0 15 6
## $ AC026464.2 : int 4 2 3 184 3 66
## $ RASA4B : int 2 0 0 59 3 0
## $ AC026954.2 : int 1 1 1 269 1 1
## $ class : Factor w/ 2 levels "EBVncp","ncp": 1 1 1 2 2 2Split into a training and testing set.
set.seed(123)
training <- ncpEBV_matrix[c(1,2,5,6),] #two each
testing <- ncpEBV_matrix[c(3,4),] #one each table(training$class)##
## EBVncp ncp
## 2 2table(testing$class)##
## EBVncp ncp
## 1 1Lets use a randomForest model to see how well it trains and predicts on unseen data.
library(randomForest)## randomForest 4.7-1.2## Type rfNews() to see new features/changes/bug fixes.str(training
)## 'data.frame': 4 obs. of 21 variables:
## $ AC139530.2 : int 520 0 1 0
## $ AC011511.4 : int 617 0 0 4
## $ AC092143.1 : int 211 0 0 0
## $ AL645465.1 : int 93 3 1 2
## $ AL353997.3 : int 1 63 1 1
## $ AC022384.1 : int 15 203 2 3
## $ MPO : int 0 36 3 0
## $ TRIM6-TRIM34: int 0 54 0 0
## $ FUT4 : int 5 4 0 1
## $ ITGA4 : int 0 13 0 1
## $ RPL23AP34 : int 0 1 1 7
## $ AC027088.3 : int 2 0 13 4
## $ MYOM3 : int 1 0 7 0
## $ CALCR : int 0 0 1 8
## $ MOXD1 : int 1 0 2 2
## $ SMAD5-AS1 : int 0 1 2 1
## $ HSPE1P7 : int 0 1 15 6
## $ AC026464.2 : int 4 2 3 66
## $ RASA4B : int 2 0 3 0
## $ AC026954.2 : int 1 1 1 1
## $ class : Factor w/ 2 levels "EBVncp","ncp": 1 1 2 2rf <- randomForest(training[,c(1:20)], training[,21],mtry=6,
importance=TRUE, confusion=T)rf$confusion## EBVncp ncp class.error
## EBVncp 2 0 0
## ncp 0 2 0The model scored 100% on training with all classes of EBV infected NCP identified and all classes of baseline NCP only identified.
Resulting output: EBVncp ncp class.error EBVncp 2 0 0 ncp 0 2 0
Now lets see how well it predicts on the hold out set or our testing set of 1 sample each.
prediction <- predict(rf,testing)result <- data.frame(predicted=prediction, actual=testing$class)
result## predicted actual
## HNE_1_MUT_LMP1_3_count EBVncp EBVncp
## HNE_1_WT_LMP1_1_count ncp ncpLooks like a very strong set of genes we selected in using a randomForest model made for classification on our data to predict with 100% accuracy the correct class of a two class model.
Resulting output:
predicted
HNE_1_WT_LMP1_1_count ncp ncp
2 rows
Thanks so much, looks like we can go ahead and add these genes to our database of pathologies. But it would be interesting to look up the actual genes used in the GOF the study used to confirm that the mutation in the H101R gene affected the immune response based on the EBVaNCP GOF.
Keep checking back.