Generating input files to use with Cellchat from OmniPathR

[saezlab]

Generating input files to use with Cellchat from OmniPathR

[deeenes]

Hi @jv-20,

Hello,

Thank you so much for this amazing resource. Great job, very useful resource +1

Thanks :)

I have few quick questions. I would like to use the underlying inter cellular network database from Omnipath to run cellchat.

A colleague already has generated CellChatDB input from OmniPath, he will probably comment and help you with this.

I am interested in interactions that are literature supported. So I applied a filter on column n_references and got a total list of 8072 unique interactions by running the following code

icn.lr_all<- import_intercell_network()
icn.lr_all$interaction_nm=paste0(icn.lr_all$source_genesymbol,":",icn.lr_all$target_genesymbol)
icn.lr_all.interaction_nm=data.frame(in_name=unique(icn.lr_all$interaction_nm))
icn.lr_all.interaction_nm$lig=sapply(strsplit(icn.lr_all.interaction_nm$in_name,split=":",fixed=T),'[',1)
icn.lr_all.interaction_nm$rec=sapply(strsplit(icn.lr_all.interaction_nm$in_name,split=":",fixed=T),'[',2)

#8072

Question 1 : Is this correct way to filter for literature supported interactions?

That's not bad, although very recently (this week) we added new features to OmnipathR to facilitate the quality filtering of intercellular networks. You can read the details in the manual, especially at import_intercell_network and filter_intercell_network: https://static.omnipathdb.org/omnipathr_manual.pdf

A basic automatic filtering:

icn_hc <- import_intercell_network(high_confidence = TRUE)

For custom parameters:

icn <- import_intercell_network()
icn_hc <- filter_intercell_network(min_references = 1, consensus_percentile = 75)

Just for convenience, we also added a parameter simplify, and a function simplify_intercell_network to get rid of some of the columns which are often not needed:

icn <- import_intercell_network()
icn_hc <- filter_intercell_network(min_references = 1, consensus_percentile = 75, simplify = TRUE)
# or:
require(magrittr)
icn_hc %<>% simplify_intercell_network

Question 2 : It would be really nice, if you could give some advice on generating input files to use with cellchat. Cellchat requires four input files. Is there function to add annotation information to a list of interactions?

In OmniPath we take all databases apart into atomic pieces of information and build our database from these. To re-create the original data structure of a database from OmniPath, you need to be careful whether the data from other sources make sense. I will comment below how these information can be recovered from OmniPath, maybe my colleague who has already done it will give you more practical comments.

* interaction_input <- CellChatDB$interaction

With this you are already on track with the import_intercell_network, as we discussed above.

* complex_input <- CellChatDB$complex

This is just a list of components for complexes, the complexes are part of the interactions in OmniPath, and their components can be easily recovered from their name (IDs separated by underscores).

* cofactor_input <- CellChatDB$cofactor

Groups of cofactors in CellChatDB are assigned to families of ligands or receptors. We process this information in a way that the cofactors are annotated in OmniPath intercell as ligand regulators or receptor regulators, while their interactions to ligands and receptors (and whether it's stimulation or inhibition) becomes part of OmniPath interactions. These information can be put together again with import_intercell_network:

cofactor_ligand <- import_intercell_network(
    transmitter_param = list(parent = 'ligand_regulator', causality = NA),
    receiver_param = list(parent = 'ligand', causality = NA)
)
cofactor_receptor <- import_intercell_network(
    transmitter_param = list(parent = 'receptor_regulator', causality = NA),
    receiver_param = list(parent = 'receptor', causality = NA)
)

* geneInfo <- CellChatDB$geneInfo

I am not sure how this is necessary for the operation of CellChatDB. Some of these gene synonyms you can recover from UniProt, using for example OmnipathR's translate_ids function.

Question 3 : In CellChat interaction_input , Ligand-Receptor interactions were assigned a pathway name. From the annotation available in omnipath, is there a column that associates interactions to pathways?

The CellChatDB pathway annotation can be recovered from OmniPath annotations, using the import_omnipath_annotations and the annotated_network functions:

cc_icn <- import_intercell_network(resources = 'CellChatDB')
cc_annot <- import_omnipath_annotations(resources = 'CellChatDB', wide = TRUE)
cc_icn_annot <- annotated_network(cc_icn, cc_annot, pathway, category)

However, you are using OmniPath because you want to use not only CellChatDB data. And to decide which other databases can be used for this purpose, depends primarily how CellChat algorithm uses this data. OmniPath has pathway annotations from a number of resources (KEGG, SignaLink, SIGNOR, NetPath, etc), all have different characteristics, and any of them can be technically used the same way as shown above.

Question 4 : Is there a function to retrieve entrez_ids for the ligands / receptors from the interaction table?

Yes, translate_ids can do this:

require(magrittr)
cc_icn <- import_intercell_network(resources = 'CellChatDB')

cc_icn %<>%
    translate_ids(
        source, source_entrez,
        uniprot, entrez,
        uploadlists = FALSE
    ) %>%
    translate_ids(
        target, target_entrez,
        uniprot, entrez,
        uploadlists = FALSE
    )

This doesn't translate complexes, for that you can use uniprot_full_id_mapping_table, and write a little function to split their names, translate and join again.

uniprot_entrez <- uniprot_full_id_mapping_table(to = 'entrez')

Thank you again for this helpful resource.

No problem, I am glad you find it useful!

[dbdimitrov]

Hi @jv-20,

Please see an issue with a similar point that I previously raised to the developers of CellChat. In case this provides you with a more detailed idea of the CellChat database and which are the relevant fields for it:
sqjin/CellChat#132

In my code, I only account for complexes and keep the original co-factors from CellChatDB. Meaning, I simply replace the interactions and complexes using those from OmniPath and reassign the to the original CellChatDB (ccDB <- CellChat::CellChatDB.human). Certainly not perfect and it could be cleaned a bit and improved. Yet, as @deeenes pointed out this information is already in OmniPath and I will also likely come back to this at some point.

Also, unfortunately, my code is specific to a certain OmniPath version and will not directly work with the information that you get from import_intercell_network.

Nevertheless, I believe the only changes that might need to be done are related to these fields (if I remember correctly):
genesymbol_intercell_source" = target, "genesymbol_intercell_target" = source, "ligand" = source_genesymbol, "receptor" = target_genesymbol

#' Helper Function to Format to CellChatDB
#' @param ccDB Inbuilt cellchatDB object
#' @param op_resource omnipath resource
#' @param exclude_anns annotations to exclude
#' @import tibble
#' @importFrom magrittr %>%
#' @importFrom stringr str_glue str_detect str_replace str_replace_all
cellchat_formatDB <- function(ccDB, op_resource, exclude_anns){
    # get complexes and interactions from omnipath
    complex_interactions <- op_resource %>%
        select(
            "genesymbol_intercell_source" = target,
            "genesymbol_intercell_target" = source,
            "ligand" = source_genesymbol,
            "receptor" = target_genesymbol,
            "evidence" = sources,
            category_intercell_source,
            category_intercell_target,
            is_directed,
            is_stimulation,
            is_inhibition
        ) %>%
        unite("annotation",
              c(category_intercell_source, category_intercell_target),
              sep="-") %>%
        unite("interaction_name", c(ligand, receptor), remove = FALSE) %>%
        mutate(pathway_name = "",
               agonist = "",
               antagonist = "",
               co_A_receptor = "",
               co_I_receptor = "") %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), ""))

    # Get OmniPath directed info
    omni_directions <- complex_interactions %>%
        select(interaction_name,
               is_directed,
               is_stimulation,
               is_inhibition,
               co_A_receptor,
               co_I_receptor
        ) %>%
        mutate(direction = pmap(., function(interaction_name,
                                            is_directed,
                                            is_stimulation,
                                            is_inhibition,
                                            co_A_receptor,
                                            co_I_receptor
        ){

            log_trace(interaction_name)

            if(co_A_receptor=="" & co_I_receptor==""){

                if(is_directed==1){
                    if(is_stimulation==1 & is_inhibition==1){
                        return("Both")
                    } else if(is_stimulation==1){
                        return("Stimulation")
                    } else if(is_inhibition==1){
                        return("Inhibition")
                    }
                }
            } else{
                return(NA)
            }
        }
        )) %>% unnest(direction)

    # get complex type of interaction from OmniPath
    omni_interactions <- complex_interactions %>%
        # set LR interactions as rowname
        left_join(., (omni_directions %>%
                          select(interaction_name,
                                 direction)),
                  by = "interaction_name") %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), "")) %>%
        mutate(co_A_receptor = ifelse(.data$co_A_receptor == "" & (direction == "Stimulation") | (direction == "both"),
                                      "Stimulation", .data$co_A_receptor),
               co_I_receptor = ifelse(.data$co_I_receptor == "" & (direction == "Inhibition") | (direction == "both"),
                                      "Inhibition", .data$co_I_receptor)) %>%
        # remove duplicates and assign to colnames
        mutate("interaction_name2" = interaction_name) %>%
        distinct_at(.vars="interaction_name2", .keep_all = TRUE) %>%
        column_to_rownames("interaction_name2") %>%
        # NOTE: ligand - (subunit_1 + subunit_2)
        mutate(interaction_name_2 = str_glue("{ligand} - {receptor}")) %>%
        mutate(interaction_name_2 = ifelse(str_detect(.data$receptor, "^COMPLEX"),
                                           str_glue("{ligand} -", "{str_split(receptor, pattern='_')}"), interaction_name_2)) %>%
        mutate(interaction_name_2 = ifelse(str_detect(.data$ligand, "^COMPLEX"),
                                           str_glue("{ligand} -", "{str_split(ligand, pattern='_')}"), interaction_name_2)) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "c", "")) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "COMPLEX:", "")) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "-", "- ")) %>%
        mutate(interaction_name_2 = str_replace_all(interaction_name_2, ", ", "+")) %>%
        mutate(interaction_name_2 = str_replace_all(interaction_name_2, '"', ""))

    # Get Omni Complexes
    omni_complexes <- complex_interactions %>%
        filter(str_detect(ligand, "_") |
                   str_detect(receptor, "_")) %>%
        select(ligand, receptor)

    # Convert to CellChat format
    omni_complexes <- union(omni_complexes$ligand,
                            omni_complexes$receptor) %>%
        enframe() %>%
        separate(col=value, sep="_",
                 into = c("subunit_1", "subunit_2",
                          "subunit_3", "subunit_4", "subunit_5"), remove=FALSE) %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), "")) %>%
        select(-name) %>%
        column_to_rownames("value")

    # Replace Default DB with OmniPath Resource
    ccDB$interaction <- omni_interactions %>%
        filter(!(annotation %in% exclude_anns))

    ccDB$complex <- omni_complexes

    return(ccDB)
}

I hope that it helps.

[dbdimitrov]

Hi Janaki,

Unfortunately, the way that I use CellChat is limited to inferring the interactions alone. As such, I did not look into the pathways, so I can't be of much help in this case.

Best wishes,
Daniel

[deeenes]

Hi Janaki,

I am glad it worked out, I am closing this issue, feel free to open another one if you have more questions.

Best,

Denes

[jv-20]

Hi @deeenes

Many thanks for your detailed answers and helping with my questions. Much appreciated. It all makes sense. I will try them out.

I look forward to your colleagues comments on how he generated CellChatDB input from OmniPath.

Thanks again, this is exciting :)

[jv-20]

Thank you so much @dbdimitrov for the reply and sharing your code. Yes, it is helpful for me. I will try it out.

Can I ask you one more question? Would it be possible to comment on how you handled the pathway column for those interactions that are not from CellChatDB?

Thank you,
Janaki

[jv-20]

Thank you Daniel

[jv-20]

Hello @deeenes

Few quick questions

First one is based on the column name "n_references" when running the function, op <- import_intercell_network() .It returns 37,647 rows and 12424 unique interactions. The minimum value of the n_references in this table is 1.

Could you please give some insights on what is the difference between column "references" and "n_references" in this table. Some interactions seem to have a value of 1 in "n_references" but no value in "references" column.

Second question is based on filter_intercell_network When I use the following code, it returns 2369 interactions.

icn_hc <- filter_intercell_network(op, min_references = 1, unique_pairs=T)

My question is , is the min_references parameter in filter_intercell_network based on "n_references" column from the results of import_intercell_network() ?

Last question is for min_curation_effort parameter. Manual says this is a minimum value of curation effort (resource-reference pairs) for network interactions. Does this mean, the interactions with min_curation_effort >=1 is manually verified?

Thank you in advance
Janaki

[deeenes]

Hi Janaki,

First one is based on the column name "n_references" when running the function, op <- import_intercell_network() .It returns 37,647 rows and 12424 unique interactions. The minimum value of the n_references in this table is 1.

Could you please give some insights on what is the difference between column "references" and "n_references" in this table. Some interactions seem to have a value of 1 in "n_references" but no value in "references" column.

Thanks for noticing this, it was a bug, counting NA as one. I fixed it in f7ed39b, will be in the repo soon.

Second question is based on filter_intercell_network When I use the following code, it returns 2369 interactions.

icn_hc <- filter_intercell_network(op, min_references = 1, unique_pairs=T)

My question is , is the min_references parameter in filter_intercell_network based on "n_references" column from the results of import_intercell_network() ?

Yes it is

Last question is for min_curation_effort parameter. Manual says this is a minimum value of curation effort (resource-reference pairs) for network interactions. Does this mean, the interactions with min_curation_effort >=1 is manually verified?

Most of the references in OmniPath are from manual curation. Some of them might be from large scale screens, and only a few resources contain literature mining. We don't have an option at the moment to perfectly separate these, but the omnipath, pathwayextra, kinaseextra and ligrecextra datasets provide an easy way to select the better curated interactions (the extra datasets are from the less curated resources).

Best,

Denes

Thank you in advance
Janaki

[jv-20]

Hi Denes,

Thank you so much for the reply, that's great.

Thanks a lot
Janaki

[superman2412]

Hi @jv-20,

Please see an issue with a similar point that I previously raised to the developers of CellChat. In case this provides you with a more detailed idea of the CellChat database and which are the relevant fields for it:
sqjin/CellChat#132

In my code, I only account for complexes and keep the original co-factors from CellChatDB. Meaning, I simply replace the interactions and complexes using those from OmniPath and reassign the to the original CellChatDB (ccDB <- CellChat::CellChatDB.human). Certainly not perfect and it could be cleaned a bit and improved. Yet, as @deeenes pointed out this information is already in OmniPath and I will also likely come back to this at some point.

Also, unfortunately, my code is specific to a certain OmniPath version and will not directly work with the information that you get from import_intercell_network.

Nevertheless, I believe the only changes that might need to be done are related to these fields (if I remember correctly):
genesymbol_intercell_source" = target, "genesymbol_intercell_target" = source, "ligand" = source_genesymbol, "receptor" = target_genesymbol

#' Helper Function to Format to CellChatDB
#' @param ccDB Inbuilt cellchatDB object
#' @param op_resource omnipath resource
#' @param exclude_anns annotations to exclude
#' @import tibble
#' @importFrom magrittr %>%
#' @importFrom stringr str_glue str_detect str_replace str_replace_all
cellchat_formatDB <- function(ccDB, op_resource, exclude_anns){
    # get complexes and interactions from omnipath
    complex_interactions <- op_resource %>%
        select(
            "genesymbol_intercell_source" = target,
            "genesymbol_intercell_target" = source,
            "ligand" = source_genesymbol,
            "receptor" = target_genesymbol,
            "evidence" = sources,
            category_intercell_source,
            category_intercell_target,
            is_directed,
            is_stimulation,
            is_inhibition
        ) %>%
        unite("annotation",
              c(category_intercell_source, category_intercell_target),
              sep="-") %>%
        unite("interaction_name", c(ligand, receptor), remove = FALSE) %>%
        mutate(pathway_name = "",
               agonist = "",
               antagonist = "",
               co_A_receptor = "",
               co_I_receptor = "") %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), ""))

    # Get OmniPath directed info
    omni_directions <- complex_interactions %>%
        select(interaction_name,
               is_directed,
               is_stimulation,
               is_inhibition,
               co_A_receptor,
               co_I_receptor
        ) %>%
        mutate(direction = pmap(., function(interaction_name,
                                            is_directed,
                                            is_stimulation,
                                            is_inhibition,
                                            co_A_receptor,
                                            co_I_receptor
        ){

            log_trace(interaction_name)

            if(co_A_receptor=="" & co_I_receptor==""){

                if(is_directed==1){
                    if(is_stimulation==1 & is_inhibition==1){
                        return("Both")
                    } else if(is_stimulation==1){
                        return("Stimulation")
                    } else if(is_inhibition==1){
                        return("Inhibition")
                    }
                }
            } else{
                return(NA)
            }
        }
        )) %>% unnest(direction)

    # get complex type of interaction from OmniPath
    omni_interactions <- complex_interactions %>%
        # set LR interactions as rowname
        left_join(., (omni_directions %>%
                          select(interaction_name,
                                 direction)),
                  by = "interaction_name") %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), "")) %>%
        mutate(co_A_receptor = ifelse(.data$co_A_receptor == "" & (direction == "Stimulation") | (direction == "both"),
                                      "Stimulation", .data$co_A_receptor),
               co_I_receptor = ifelse(.data$co_I_receptor == "" & (direction == "Inhibition") | (direction == "both"),
                                      "Inhibition", .data$co_I_receptor)) %>%
        # remove duplicates and assign to colnames
        mutate("interaction_name2" = interaction_name) %>%
        distinct_at(.vars="interaction_name2", .keep_all = TRUE) %>%
        column_to_rownames("interaction_name2") %>%
        # NOTE: ligand - (subunit_1 + subunit_2)
        mutate(interaction_name_2 = str_glue("{ligand} - {receptor}")) %>%
        mutate(interaction_name_2 = ifelse(str_detect(.data$receptor, "^COMPLEX"),
                                           str_glue("{ligand} -", "{str_split(receptor, pattern='_')}"), interaction_name_2)) %>%
        mutate(interaction_name_2 = ifelse(str_detect(.data$ligand, "^COMPLEX"),
                                           str_glue("{ligand} -", "{str_split(ligand, pattern='_')}"), interaction_name_2)) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "c", "")) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "COMPLEX:", "")) %>%
        mutate(interaction_name_2 = str_replace(interaction_name_2, "-", "- ")) %>%
        mutate(interaction_name_2 = str_replace_all(interaction_name_2, ", ", "+")) %>%
        mutate(interaction_name_2 = str_replace_all(interaction_name_2, '"', ""))

    # Get Omni Complexes
    omni_complexes <- complex_interactions %>%
        filter(str_detect(ligand, "_") |
                   str_detect(receptor, "_")) %>%
        select(ligand, receptor)

    # Convert to CellChat format
    omni_complexes <- union(omni_complexes$ligand,
                            omni_complexes$receptor) %>%
        enframe() %>%
        separate(col=value, sep="_",
                 into = c("subunit_1", "subunit_2",
                          "subunit_3", "subunit_4", "subunit_5"), remove=FALSE) %>%
        mutate_at(vars(everything()), ~ replace(., is.na(.), "")) %>%
        select(-name) %>%
        column_to_rownames("value")

    # Replace Default DB with OmniPath Resource
    ccDB$interaction <- omni_interactions %>%
        filter(!(annotation %in% exclude_anns))

    ccDB$complex <- omni_complexes

    return(ccDB)
}

I hope that it helps.

Hi, thank you so much for writing this function! Unfortunately, it just doesn't seem to work. The output after running this function is exactly the same as the original CellChatDB. Probably because the function is not compatible with import_intercell_network?

I extracted all interactions from Omnipath by:

all_interactions <- import_intercell_network ()

And saved it as a .csv file (attached) below:

all_interactions.csv

I would be really grateful to you if you could please help me add these interactions to CellChatDB?

Also, CellChat authors have published this tutorial as to adding interactions here Probably this will be of help?

[dbdimitrov]

Hi,

You can use liana - our recent package to do this, with the following code:

require(liana)
# Get custom OmniPath Resource
cust_omni <- generate_omni(loc_consensus_percentile = 30,
                           consensus_percentile = NULL,
                           transmitter_topology = c('secreted',
                                                    'plasma_membrane_transmembrane',
                                                    'plasma_membrane_peripheral'),
                           receiver_topology = c('plasma_membrane_transmembrane',
                                                 'plasma_membrane_peripheral'),
                           min_curation_effort = 1,
                           ligrecextra = FALSE,
                           remove_complexes = FALSE,
                           simplify = FALSE
)
cust_omni %>% glimpse

# test
cellchatDB_Omni <- liana:::cellchat_formatDB(ccDB = CellChat::CellChatDB.human,
                                             op_resource = ligrec,
                                             exclude_anns = c())
cellchatDB_Omni %>% glimpse

This vignette provides further detail about whats under the hood - i.e. OmniPathR : https://saezlab.github.io/liana/articles/liana_custom_op.html

I hope this helps!

Regards,
Daniel

[dbdimitrov]

Perhaps, it would be even easier to follow the same tutorial to simply generate a resource according to your preference and then use liana to call Cellchat with it

[superman2412]

@dbdimitrov Thank you so much for your prompt reply! It works! I have 2 quick questions:

1. Does using this code add to the existing CellChatDB or does it replace the contents of the CellChatDB? If it replaces the contents of CellChatDB, how should I customise the code so that LIANA can add contents of Omnipath on top of what is already there in CellChatDB.human?

2. Ideally, it would be great to have as much coverage of the interactions as much possible. Do you think that for parameter opresource resource <- import_intercell_network() would lead the greatest number of possible interactions?

Really grateful to you for your help in this.

[dbdimitrov]

Hi,

1. It replaces them, but I would presume that large parts of CellChatDB are in the OmniPath CCC resource.
2. I'm not too certain if this would work, perhaps modifying the filtering parameters from the vignette to be as non-selective as possible would be the easiest way to get the highest coverage.

[superman2412]

Hi,

1. It replaces them, but I would presume that large parts of CellChatDB are in the OmniPath CCC resource.
2. I'm not too certain if this would work, perhaps modifying the filtering parameters from the vignette to be as non-selective as possible would be the easiest way to get the highest coverage.

Got it thank you so much! Actually using resource <- import_intercell_network() worked till the steps shown in the three screenshots below.

However, as below there's a problem when I try to run, cellchat <- computeCommunProbPathway(object = cellchat). Any suggestions as to how I might be able to fix that?

[superman2412]

Hi @dbdimitrov , I also tried running cellchat with the code you sent.

require(liana)
# Get custom OmniPath Resource
cust_omni <- generate_omni(loc_consensus_percentile = 30,
                           consensus_percentile = NULL,
                           transmitter_topology = c('secreted',
                                                    'plasma_membrane_transmembrane',
                                                    'plasma_membrane_peripheral'),
                           receiver_topology = c('plasma_membrane_transmembrane',
                                                 'plasma_membrane_peripheral'),
                           min_curation_effort = 1,
                           ligrecextra = FALSE,
                           remove_complexes = FALSE,
                           simplify = FALSE
)
cust_omni %>% glimpse

# test
cellchatDB_Omni <- liana:::cellchat_formatDB(ccDB = CellChat::CellChatDB.human,
                                             op_resource = ligrec,
                                             exclude_anns = c())
cellchatDB_Omni %>% glimpse

I only made a minor adjustment I changed op_resource = ligrec into op_resource = cust_omni .

And its the same story, it ran till the first three screenshots that I shared (as above) but it didn't run the code, with exactly the same error (see the screenshot below). I was wondering if we could do something to fix this?

Really appreciate your help! Thank you so much!

[superman2412]

@dbdimitrov Would be really grateful to you if you could shed some light on this issue.

[dbdimitrov]

Hi @superman2412,

The error that you get is unfortunately expected, since we don't modify the pathway information in the CellChatDB.

OmniPathR contains the information that you might need to modify the DB appropriately under import_annotations. However, I don't see an easy way to do that and have not tried it myself. So, I can be of little help.

I hope this helps.
Regards,
Daniel