Body: <p>I am willing to reproduce part of the analysis from "The human transcriptome across tissues and individuals" (Melé et all, 2015).</p>
  
  <p>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). </p>
  
  <p>In my process I did:</p>
  
  <ol>
  <li>Load the FPKC matrix</li>
  <li>Remove those transcrits not having a level of 0.1 in minim 10 samples</li>
  <li>log10 transform after adding 1 to the matrix (to avoid 0s transformation) </li>
  <li>Quantile normalize the matrix</li>
  <li>Compute the distance matrix between samples</li>
  <li>MDS and plot</li>
  </ol>
  
  <p>But my plot has nothing to to with the one in the paper. </p>
  
  <p><a href="https://i.sstatic.net/gyXCl.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/gyXCl.png" alt="MDS from GTEx v6 transcriptome data"></a></p>
  
  <p>Can anyone see any error in my process o see what I am doing that is totally different from GTEx process pipeline?</p>
  
  <hr>
  
  <h2>Extra info</h2>
  
  <pre><code>library( preprocessCore )
  library( smacof )
  
  # Load gene expression matrix from file
  gen &lt;- 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 &lt;- gen[ , c(1,2) ]
  gen_clean &lt;- as.matrix( gen[ , -seq( 2 ) ] )
  
  # Load map file
  sample_map &lt;- read.delim( "data/v6/GTEx_Data_V6_Annotations_SampleAttributesDS.txt" )
  sample_map$SAMPID &lt;- gsub( "-", "\\.", sample_map$SAMPID )
  rownames( sample_map ) &lt;- sample_map$SAMPID
  
  # Match data-sets
  common_samples &lt;- intersect( 
      colnames( gen_clean ), rownames( sample_map )
  )
  gen_clean &lt;- gen_clean[ , common_samples ]
  sample_map &lt;- sample_map[ common_samples, ]
  
  # Select tissues
  tissue_to_get &lt;- c( "muscle", "brain", "heart", "blood", "pituitary", "adipose tissue", "nerve", "lung", "breast", "blood vessel", "thyroid", "liver", "skin", "testis" )
  
  samples_to_use &lt;- sample_map$SAMPID[ 
       tolower( as.character( sample_map$SMTS ) ) %in% tissue_to_get 
  ]
  
  # Filter gene expression matrix according to GTEx
  trs_to_use &lt;- apply( X = gen_clean, MARGIN = 1, FUN = function( row, exp_th, min_sam ) {
    sum( row &gt; exp_th ) &gt;= min_sam
  }, exp_th = 0.1, min_sam = 10 )
  gen_rd &lt;- gen_clean[ trs_to_use, samples_to_use  ]
  
  # log10 transformation and normalization
  gen_norm &lt;- normalize.quantiles( log10( gen_rd[ , list_samples ] + 1 ) )
  
  # MDS
  d &lt;- dist( t( gen_norm ) )
  mds &lt;- mds( d, type = "ratio", ndim = 3 )
  [ ... creation of color map ... ]
  plot( x = mds$conf[ , 1 ], y = mds$conf[ , 3 ], col = clr2, pch = 20 )
  </code></pre>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-10T11:24:37.473Z
Comments.Text: 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
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-10T13:18:30.207Z
Comments.Text: 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.
Comments.UserDisplayName: carlesh
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2018-08-10T14:05:21.923Z
Comments.Text: 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.
Comments.UserDisplayName: llrs
Tags: <r><rna-seq><transcriptome><figure-reproduction>
Title: Reproducing GTEx transcriptome analysis
ThreadText: # Reproducing GTEx transcriptome analysis
  
  (Asked by: username_0 on 2018-08-08)
  
  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](https://i.sstatic.net/gyXCl.png)](https://i.sstatic.net/gyXCl.png)
  
  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 )
  
  ```
  
  **Tags:** r, rna-seq, transcriptome, figure-reproduction
  
  > **Comments:**
  > **username_1**: 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  
  > **username_0**: 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.  
  > **username_1**: 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.  
  
  ---
  
  ## 0 answers
  
  