Skip to content

Commit eae1f4c

Browse files
committed
pairwise comparison and description for genes
1 parent 1351500 commit eae1f4c

File tree

1 file changed

+15
-16
lines changed

1 file changed

+15
-16
lines changed

R/Analyze_genecatalog.Rmd

+15-16
Original file line numberDiff line numberDiff line change
@@ -553,10 +553,19 @@ library(DESeq2)
553553

554554
Prepare Data:
555555

556+
```{r}
557+
Description <- annotations %>%
558+
select(ko_id, kegg_hit) %>%
559+
distinct(ko_id,.keep_all = TRUE) %>%
560+
column_to_rownames("ko_id") %>%
561+
as.data.frame()
562+
````
563+
556564
```{r }
557565
# Assuming annotation_cpm is your count data (CPM)
558566
dds <- DESeqDataSetFromMatrix(countData = annotation_genecounts,
559567
colData = metadata[colnames(annotation_genecounts),],
568+
rowData = Description,
560569
design = ~ group)
561570
```
562571

@@ -569,9 +578,12 @@ Normalize the data using the DESeq2 normalization methods, and estimate differen
569578
```
570579

571580
```{r}
572-
res <- results(dds) %>% as.data.frame()
581+
res <- results(dds,contrast=c(group_variable,"AD","HC")) %>% as.data.frame()
582+
583+
res %>% merge( Description, by = "row.names", all = TRUE) %>%
584+
column_to_rownames("Row.names") -> res
573585
574-
head(res %>% arrange(padj),n=10)
586+
res %>% arrange(padj) %>% select(all_of(c("baseMean", "log2FoldChange", "padj","kegg_hit"))) %>% head(n=10)
575587
```
576588

577589
```{r }
@@ -588,6 +600,7 @@ head(res %>% arrange(padj),n=10)
588600
```
589601

590602

603+
591604
## Alternative ways to normalize
592605

593606
An other interesting way to normalize would be to normalize to single copy genes.
@@ -596,21 +609,7 @@ Example the one from [MuSiCC](https://github.com/borenstein-lab/MUSiCC/blob/mast
596609
MUSiCC: A marker genes based framework for metagenomic normalization and accurate profiling of gene abundances in the microbiome. Ohad Manor and Elhanan Borenstein. Genome Biology
597610

598611

599-
## Integrate them in a SummarizedExperiment object
600-
601-
```{r }
602612

603-
SE <- SummarizedExperiment(
604-
assays = list(gcpm = annoation_cpm),
605-
colData = metadata[colnames(annoation_cpm), ],
606-
#rowData = annotations
607-
)
608-
609-
```
610-
611-
```{}
612-
library(tidybulk)
613-
```
614613

615614

616615

0 commit comments

Comments
 (0)