Answers.Body: <p>I'd use <code>dplyr</code> for this.</p>
  
  <p>First, <code>mouse paralogues</code> contains no values for <code>mmusculus_paralog_associated_gene_name</code> that are NA:</p>
  
  <pre><code>library(dplyr)
  mouse_paralogues %&gt;% 
    filter(is.na(mmusculus_paralog_associated_gene_name)) %&gt;% 
    nrow()
  
  [1] 0
  </code></pre>
  
  <p>But there are values which are empty strings, so you can filter for those if you wish.</p>
  
  <pre><code>count_mouse_paralogues &lt;- mouse_paralogues %&gt;%
    filter(hsapiens_homolog_associated_gene_name != "") %&gt;%
    count(external_gene_name)
  </code></pre>
  
  <p>Confirm that this gives the expected result (n = 4) for Plin4:</p>
  
  <pre><code>count_mouse_paralogues %&gt;%
    filter(external_gene_name == "Plin4")
  
  # A tibble: 1 x 2
    hsapiens_homolog_associated_gene_name     n
    &lt;chr&gt;                                 &lt;int&gt;
  1 Plin4                                     4
  </code></pre>
  
  <p>And check <code>mouse_paralogues</code>:</p>
  
  <pre><code>mouse_paralogues %&gt;% 
    filter(hsapiens_homolog_associated_gene_name == "PLIN4")
  
    hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
  1                                 PLIN4              Plin4                                  Plin3
  2                                 PLIN4              Plin4                                  Plin5
  3                                 PLIN4              Plin4                                  Plin2
  4                                 PLIN4              Plin4                                  Plin1
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-31T05:57:17.057Z
Answers.LastActivityDate: 2018-01-31T06:04:45.530Z
Answers.LastEditDate: 2018-01-31T06:04:45.530Z
Answers.UserDisplayName: neilfws
Body: <p>I am trying to count the number of paralogues for the mouse homologues of the human protein-coding genes using BioMart. But for example in the 'PLIN4' gene its counting 35,000 paralogues instead of 4. </p>
  
  <p>We think it is because some genes have one to many paralogues which causes repeats. When I run a single gene its gives me back the correct number of paralogues. Is there a way to either remove these repeats from the results or a way around this so that BioMart doesn't output these repeats. </p>
  
  <p>I have also thought of maybe running one gene at a time, then counting it by setting up some sort of loop so that it does all of the genes from the list automatically. </p>
  
  <p>The code I have written so far is:</p>
  
  <pre><code># Load the biomaRt package:
  
  library(biomaRt)
  ensembl_hsapiens &lt;- useMart("ensembl", 
                            dataset = "hsapiens_gene_ensembl")
  ensembl_mouse &lt;- useMart("ensembl", 
                         dataset = "mmusculus_gene_ensembl")
  
  # Get all human protein coding genes:
  
  hsapien_PC_genes &lt;- getBM(attributes = c("ensembl_gene_id", 
                                           "external_gene_name"), 
                            filters = "biotype", 
                            values = "protein_coding", 
                            mart = ensembl_hsapiens)
  
  
  ensembl_gene_ID &lt;- hsapien_PC_genes$ensembl_gene_id
  
  # Get mouse homologues
  
  mouse_homologues &lt;- getBM(attributes = c("ensembl_gene_id", "external_gene_name", 
                                         "mmusculus_homolog_associated_gene_name"), 
                          filters = "ensembl_gene_id", 
                          values = c(ensembl_gene_ID), 
                          mart = ensembl_hsapiens)
  
  # Get mouse external gene name 
  
  mouse_homologues_external_gene_names &lt;- mouse_homologues$mmusculus_homolog_associated_gene_name
  
  
  mouse_paralogues &lt;- getBM(attributes = c("hsapiens_homolog_associated_gene_name",
                                         "external_gene_name",
                                         "mmusculus_paralog_associated_gene_name"), 
                          filters = "external_gene_name", 
                          values = c(mouse_homologues_external_gene_names) , mart = ensembl_mouse)
  
  # Remove genes with no paralogues 
  mouse_paralogs_data &lt;- mouse_paralogues[!(is.na(mouse_paralogues$mmusculus_paralog_associated_gene_name)
                                            | 
  mouse_paralogues$mmusculus_paralog_associated_gene_name==""), ]
  
  # Count paralogues per gene
  
  library(plyr)
  count_mouse_paralogues &lt;- count(mouse_paralogs_data, "external_gene_name")
  View(count_mouse_paralogues)
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-01-31T22:01:15.037Z
Comments.Text: Related post at SO: https://stackoverflow.com/questions/48468052/counting-the-number-of-paralogues-for-mouse-genes-gives-me-the-wrong-frequency
Comments.UserDisplayName: zx8754
Tags: <r><bioconductor><biomart><ensembl>
Title: Counting the number of paralogues for mouse genes gives me the wrong frequency in R
ThreadText: # Counting the number of paralogues for mouse genes gives me the wrong frequency in R
  
  (Asked by: username_0 on 2018-01-27)
  
  I am trying to count the number of paralogues for the mouse homologues of the human protein-coding genes using BioMart. But for example in the 'PLIN4' gene its counting 35,000 paralogues instead of 4.
  
  We think it is because some genes have one to many paralogues which causes repeats. When I run a single gene its gives me back the correct number of paralogues. Is there a way to either remove these repeats from the results or a way around this so that BioMart doesn't output these repeats.
  
  I have also thought of maybe running one gene at a time, then counting it by setting up some sort of loop so that it does all of the genes from the list automatically.
  
  The code I have written so far is:
  
  ```
  # Load the biomaRt package:
  
  library(biomaRt)
  ensembl_hsapiens <- useMart("ensembl", 
                            dataset = "hsapiens_gene_ensembl")
  ensembl_mouse <- useMart("ensembl", 
                         dataset = "mmusculus_gene_ensembl")
  
  # Get all human protein coding genes:
  
  hsapien_PC_genes <- getBM(attributes = c("ensembl_gene_id", 
                                           "external_gene_name"), 
                            filters = "biotype", 
                            values = "protein_coding", 
                            mart = ensembl_hsapiens)
  
  ensembl_gene_ID <- hsapien_PC_genes$ensembl_gene_id
  
  # Get mouse homologues
  
  mouse_homologues <- getBM(attributes = c("ensembl_gene_id", "external_gene_name", 
                                         "mmusculus_homolog_associated_gene_name"), 
                          filters = "ensembl_gene_id", 
                          values = c(ensembl_gene_ID), 
                          mart = ensembl_hsapiens)
  
  # Get mouse external gene name 
  
  mouse_homologues_external_gene_names <- mouse_homologues$mmusculus_homolog_associated_gene_name
  
  mouse_paralogues <- getBM(attributes = c("hsapiens_homolog_associated_gene_name",
                                         "external_gene_name",
                                         "mmusculus_paralog_associated_gene_name"), 
                          filters = "external_gene_name", 
                          values = c(mouse_homologues_external_gene_names) , mart = ensembl_mouse)
  
  # Remove genes with no paralogues 
  mouse_paralogs_data <- mouse_paralogues[!(is.na(mouse_paralogues$mmusculus_paralog_associated_gene_name)
                                            | 
  mouse_paralogues$mmusculus_paralog_associated_gene_name==""), ]
  
  # Count paralogues per gene
  
  library(plyr)
  count_mouse_paralogues <- count(mouse_paralogs_data, "external_gene_name")
  View(count_mouse_paralogues)
  
  ```
  
  **Tags:** r, bioconductor, biomart, ensembl
  
  > **Comments:**
  > **username_1**: Related post at SO: https://stackoverflow.com/questions/48468052/counting-the-number-of-paralogues-for-mouse-genes-gives-me-the-wrong-frequency  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_2 on 2018-01-31. Score: 1
  
  I'd use `dplyr` for this.
  
  First, `mouse paralogues` contains no values for `mmusculus_paralog_associated_gene_name` that are NA:
  
  ```
  library(dplyr)
  mouse_paralogues %>% 
    filter(is.na(mmusculus_paralog_associated_gene_name)) %>% 
    nrow()
  
  [1] 0
  
  ```
  
  But there are values which are empty strings, so you can filter for those if you wish.
  
  ```
  count_mouse_paralogues <- mouse_paralogues %>%
    filter(hsapiens_homolog_associated_gene_name != "") %>%
    count(external_gene_name)
  
  ```
  
  Confirm that this gives the expected result (n = 4) for Plin4:
  
  ```
  count_mouse_paralogues %>%
    filter(external_gene_name == "Plin4")
  
  # A tibble: 1 x 2
    hsapiens_homolog_associated_gene_name     n
    <chr>                                 <int>
  1 Plin4                                     4
  
  ```
  
  And check `mouse_paralogues`:
  
  ```
  mouse_paralogues %>% 
    filter(hsapiens_homolog_associated_gene_name == "PLIN4")
  
    hsapiens_homolog_associated_gene_name external_gene_name mmusculus_paralog_associated_gene_name
  1                                 PLIN4              Plin4                                  Plin3
  2                                 PLIN4              Plin4                                  Plin5
  3                                 PLIN4              Plin4                                  Plin2
  4                                 PLIN4              Plin4                                  Plin1
  
  ```
  
  
  ---
  
  