2
$\begingroup$

I am willing to reproduce part of the analysis from "The human transcriptome across tissues and individuals" (Melé et all, 2015).

I downloaded GTEx v6 FPKM data in txt format from GTEx portal. I want to mimic the MDS plot from Supplementary Figure 3 (or get the equivalent PCA).

In my process I did:

  1. Load the FPKC matrix
  2. Remove those transcrits not having a level of 0.1 in minim 10 samples
  3. log10 transform after adding 1 to the matrix (to avoid 0s transformation)
  4. Quantile normalize the matrix
  5. Compute the distance matrix between samples
  6. MDS and plot

But my plot has nothing to to with the one in the paper.

MDS from GTEx v6 transcriptome data

Can anyone see any error in my process o see what I am doing that is totally different from GTEx process pipeline?


Extra info

library( preprocessCore )
library( smacof )

# Load gene expression matrix from file
gen <- read.delim( "data/v6/GTEx_Analysis_v6p_RNA-seq_RNA-SeQCv1.1.8_gene_rpkm.gct", skip = 2, heade = TRUE )

## First two columns corresponds to "Transcript name" and "Gene name"
gen_ann <- gen[ , c(1,2) ]
gen_clean <- as.matrix( gen[ , -seq( 2 ) ] )

# Load map file
sample_map <- read.delim( "data/v6/GTEx_Data_V6_Annotations_SampleAttributesDS.txt" )
sample_map$SAMPID <- gsub( "-", "\\.", sample_map$SAMPID )
rownames( sample_map ) <- sample_map$SAMPID

# Match data-sets
common_samples <- intersect( 
    colnames( gen_clean ), rownames( sample_map )
)
gen_clean <- gen_clean[ , common_samples ]
sample_map <- sample_map[ common_samples, ]

# Select tissues
tissue_to_get <- c( "muscle", "brain", "heart", "blood", "pituitary", "adipose tissue", "nerve", "lung", "breast", "blood vessel", "thyroid", "liver", "skin", "testis" )

samples_to_use <- sample_map$SAMPID[ 
	    tolower( as.character( sample_map$SMTS ) ) %in% tissue_to_get 
]

# Filter gene expression matrix according to GTEx
trs_to_use <- apply( X = gen_clean, MARGIN = 1, FUN = function( row, exp_th, min_sam ) {
  sum( row > exp_th ) >= min_sam
}, exp_th = 0.1, min_sam = 10 )
gen_rd <- gen_clean[ trs_to_use, samples_to_use  ]

# log10 transformation and normalization
gen_norm <- normalize.quantiles( log10( gen_rd[ , list_samples ] + 1 ) )

# MDS
d <- dist( t( gen_norm ) )
mds <- mds( d, type = "ratio", ndim = 3 )
[ ... creation of color map ... ]
plot( x = mds$conf[ , 1 ], y = mds$conf[ , 3 ], col = clr2, pch = 20 )
$\endgroup$
3
  • $\begingroup$ Welcome to the site Carles! Are you comparing your figure with the figure 3 A) or Figure 1 B) (left)? The step 4 you mention is not used in the article, also they filter without any consideration of the number of samples it is present $\endgroup$
    – llrs
    Aug 10 '18 at 11:24
  • $\begingroup$ Hi @Llopis, thanks for the welcome! I am comparing my MDS with the one from Fig S3 (Multi dimensional scaling (MDS) of GTEx samples based on gene expression). I agree that step 4 (quantile normalization) is not described but I assume they apply some kind of normalization. And they apply quantile normalization in the array set, that is the reason I applied it. $\endgroup$
    – carlesh
    Aug 10 '18 at 13:18
  • $\begingroup$ Ok, then my other suggestion is to do the distance metric as they do, 1 - Pearson correlation as they did, maybe you get closer results. $\endgroup$
    – llrs
    Aug 10 '18 at 14:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.