How to identify tissue-specific genes

identify_tissue-specific_genes_using_median_expression.r

# Aim: Identifying genes that are at least 2 fold higher mRNA levels (FPKM) in a particular tissue as 
# -compared to all other tissues. 

# Author: Gireesh Bogu
# Date: Jun 25th, 2017
# Location: CRG, Barcelona

# Problem: How to idenitfy tissue-specific genes especially when you have 
# - large number of tissues (>50) and even larger number of samples per tissue (>100 for example).
# GTEx (2) has 53 tissue sites and each tissue site has 10 to 400 samples

# Solution: Comparing mean expression leads to biased results. Here we use a median based approach where we calcualte median
# - expression of each gene per tissue and then compare the tissue with maximum median expression against the tissue that 
# - has second highest median exppression. If the difference between them is more than 2 fold then the gene will be 
# - selected as tissue-specific. (Using >=2 fold change of median expression) 

# pros: Very simple and very fast (8000 columns and 870 rows data matrix took less than a minute)
# cons: This script only finds tissue-specific or in other words tissue-enriched but not other categories as listed below.
#   Category definition like this improves the understanding of tissue-specificity even better (3)
#   1. Tissue enriched = At least fivefold higher mRNA levels (FPKM) in a particular tissue as compared to all other tissues
#   2. Group enriched = At least fivefold higher mRNA levels in a group of tissues 
#   3. Enhanced = At least fivefold higher mRNA levels in a particular tissue as compared to the average levels in all tissues
#   4. House keeping = Expressed in all	Detected in all tissues	
#   5. Mixed	= Detected in at least two tissues, but not in all, and not part of any of the categories above	
#   6. Repressed = Not detected	Not present in any of the analyzed tissues

# No statistical association has been performed as they show extremely significant pvalues (ex: ANOVA) because of large sample sizes.
# testing "Kruskal-Wallis Test" (rank based)

# important:
# 1. make sure you rename the tissue sites into major tissues (all brain types into one). : if you do not you might not detect the ones that are specific to brain tissue and have similar levels of expression in at least two tissue sites 
# 2. You can filter out the tissue-specific genes from the output.txt further by usign minimum expression threhold > 1 (to reduce noise). 
# 3. Use 5 fold for protein-coding and 2 fold for lincRNA or transposable elements. 


# input.txt
aa_srrids aa_tissues  gene1  gene2  gene3  genen
srrx1 brain 2 8 6 5
srrx2 brain 5 7 9 9
srrx3 brain 6 7 8 9
srrx4 brain 8 8 8 8
srrx5 heart 7 7 7 7
srrx6 heart 10  10  10  10
srrx7 lung  3 12  1 48
srrx8 lung  5 5 24  48
srrx9 lung  12  3 12  48
srrx10  lung  2 3 3 44

# output.txt
lung    genen    48

# code

library(dplyr)
library(tidyr)

a <- read.table("input.txt", head=T)

b <- a %>%
  group_by(aa_tissues) %>% #gropus by tissue names
  summarise_each(funs(median(., na.rm=TRUE)), -aa_srrids) #'-' ignores specific column from the data

c<- b %>% # transpose the tibble so that genes will be in rows and tissues will be in columns
  gather(var, val, 2:ncol(b)) %>%
  spread(aa_tissues, val)

d <- c %>% # Identify genes with max tissue median expression
  gather(k, v, -var) %>%
  arrange(-v, var) %>%
  group_by(var) %>%
  mutate(logi = v / lead(v, default = min(v)) >= 2) %>% # fold change >=2 (2 times greater than the next max)
  filter(v == first(v) & logi) %>%
  select(-logi)

write.table(d, file="output.txt", quote=F, sep="\t")

# helpful links:
# (1) http://www.proteinatlas.org/humanproteome/tissue+specific
# (2) https://www.gtexportal.org/home/
# (3) https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4848759/
# http://tidyverse.org/
# https://stackoverflow.com/questions/25198442/how-to-calculate-mean-median-per-group-in-a-dataframe-in-r
# https://stackoverflow.com/questions/21982987/mean-per-group-in-a-data-frame
# https://stackoverflow.com/questions/9723208/aggregate-summarize-multiple-variables-per-group-i-e-sum-mean-etc
# https://stackoverflow.com/questions/44745756/how-to-identify-the-maximum-value-that-is-2-times-higher-than-next-maximum-value#44745936