PubMA Exercise.3

From BITS wiki
Jump to: navigation, search

Clustering using the GEO Dataset browser (only for data with attached GDS ID)

GDSbrowser_logo.png

[ Main_Page | Hands-on Analysis of public microarray datasets | PubMA_Exercise.2 | PubMA_Exercise.3 |
| PubMA_Exercise.4 ]


Introduction

The GEO Dataset Browser [1] takes care of providing users clustered data and heatmaps derived from manual curation of some of the GEO datasets selected by the NCBI team. This processing is otherwise laborious and requires [R] skills that not every scientist can build. More info about this toolset can be found on the NCBI help page[2]

The web interface allows finding co-expressed genes in clusters. Co-expressed genes often belong to common signaling pathways or result from transcriptional co-regulation by common transcription factors.

The GEO Dataset browser

The remaining of this page shows key view obtained by navigating in the GEO Dataset browser Web interface

Start the tool

Start by searching NCBI's Bioproject or GEO Datasets database using the ID GSE6943. On the resulting Bioproject page you can follow the link to the GEO data sets that correspond to this microarray experiment.


Bioproject.png


Select GEO Dataset Browser data

The above link brings you to the GEO DataSets page shown next. The first reference links to GEO2R to identify DE genes while the second refers to the GEO Dataset Browser for clustering and creation of a heat map.


02-GEO2R-GDSB.png


The Analyze DataSet link opens a new tab


03-GDS-browser.png


You can find a full description of GEO DataSet records on NCBI's About GEO Data Sets page.

Download the full data table

The Download section in the right menu allows you to download the full normalized data set (Series Matrix file) in various formats

03b_Download_data.png

If you do the download, the resulting text file contains a 49 lines header which may pose problems in Excel

See how to split the file in two using some CLI magic

# download the file from the 'DataSet SOFT file' link in the main window
wget ftp://ftp.ncbi.nlm.nih.gov/geo/datasets/GDS3nnn/GDS3224/soft/GDS3224.soft.gz
 
# decompress and replace line ends bby valid unix CR
gunzip -c GDS3224.soft.gz | tr "\r" "\n" > GDS3224.soft.txt
 
# find table line#1
grep -n 'dataset_table_begin' GDS3224.soft.txt
49:!dataset_table_begin
 
# we need only the number 49!
header_end=$(grep -n 'dataset_table_begin' GDS3224.soft.txt | awk 'BEGIN{FS=":"}{print $1}')
 
# show the result
echo ${header_end}
49
 
# split the file in two
head -${header_end} GDS3224.soft.txt > GDS3224.data-header.txt
cat GDS3224.soft.txt | sed -e "1,${header_end}d" > GDS3224.data.txt
 
# inspect results
grep ^#GSM GDS3224.data-header.txt
#GSM160089 = Value for GSM160089: Diaphragm 1; src: Diaphragm
#GSM160090 = Value for GSM160090: Diaphragm 2; src: Diaphragm
#GSM160091 = Value for GSM160091: Diaphragm 3; src: Diaphragm
#GSM160092 = Value for GSM160092: Diaphragm 4; src: Diaphragm
#GSM160093 = Value for GSM160093: Diaphragm 5; src: Diaphragm
#GSM160094 = Value for GSM160094: Diaphragm 6; src: Diaphragm
#GSM160095 = Value for GSM160095: Heart 1; src: Heart (left vent)
#GSM160096 = Value for GSM160096: Heart 2; src: Heart (left vent)
#GSM160097 = Value for GSM160097: Heart 3; src: Heart (left vent)
#GSM160098 = Value for GSM160098: Heart 4; src: Heart (left vent)
#GSM160099 = Value for GSM160099: Heart 5; src: Heart (left vent)
#GSM160100 = Value for GSM160100: Heart 6; src: Heart (left vent)
 
head GDS3224.data.txt | column -t
ID_REF      IDENTIFIER  GSM160089  GSM160090  GSM160091  GSM160092  GSM160093  GSM160094  GSM160095  GSM160096  GSM160097  GSM160098  GSM160099  GSM160100
1367452_at  Sumo2       2532.9     2518.6     2384.6     2304       2360       2482.8     3166       2938.9     2953.3     2558.8     3043.3     2711.5
1367453_at  Cdc37       3464.2     3197.4     3487.1     3133.2     3432.5     3486.9     3860.2     3429.4     3381.2     4131.3     4364.2     3605.2
1367454_at  Copb2       1620.8     1870.5     1538.6     1334       1502.9     1520.3     1849.8     1852       1858.5     1483.3     1766.1     1500.5
1367455_at  Vcp         5512.5     4103.9     5746.5     4393.6     5870.7     5851.2     5408.8     4682.6     4734.7     4403.7     3940       3674.5
1367456_at  Ube2d3      6090.8     5352.2     5614.9     5249.6     5834.6     5915.9     3995.3     4356.9     4675.1     4570.4     4994.8     4231.6
1367457_at  Becn1       1093.9     1134.3     736.4      1219       774.9      712.2      892.1      998.9      782.3      710.8      923.2      975
1367458_at  Lypla2      347.8      223.9      261.4      338.8      249.6      363.7      422.2      409.9      273.3      492.2      458        403.4
1367459_at  Arf1        7665.8     7415.9     7075.9     7349.4     6406.7     6664.6     10400.5    9729.2     9679.2     9996.8     9783.7     8333
1367460_at  Gdi2        3155.7     2946.9     3589.7     3487.4     3131.5     3338       3198       2991.3     2781       2754.1     2406.9     2609.5


Walking through the tools

Several built-in tools are ready to serve you on demand. We provide below a rapid overview of the results obtained by each tool.

Find Genes

Allows users interested in a particular gene to get its expression values across samples.

03c_Find-genes.png

Compare two sets of samples

This allows selecting two groups of samples that are compared using a t-test.


04-compare2.png

You can choose between:

  • a two-tailed comparison: you check if the two groups have different expression
  • a one-tailed comparison: you check if expression in group 1 is higher than in group 2
  • a one-tailed comparison: you check if expression in group 1 is smaller than in group 2
  • value means difference test: you check if the mean fold change is exceeding a threshold, e.g. more than two-fold difference between the groups
  • rank means difference test: test based on ranks of fold changes, you test if the mean rank difference is exceeding a threshold, e.g. more than two-fold

Remember that ordinary t-tests are not well-suited for microarray data because the assumptions of the t-test are not met (equal variances between the groups), microarray experiments are very noisy and there are too few biological replicates to make accurate estimates of means and variations.

Ideally, you use a modified t-test (as in limma) to compare expression levels. In this tool of the GEO browser, they use a regular t-test without multiple testing correction (as far as I can conclude from the very concise documentation). You can set the significance level and not a FDR. This means that if you set the significance level to the lowest value (0,01) you allow 1% false positives. The rat array represents around 10000 genes so you'll do 10000 t-tests. With an alpha level of 0,01 this will generate 100 false positives.

Although the t-test in this tool is far from ideal, the two other tests also have their limitations, mainly the fact that they rely on arbitrary thresholds (3-fold? 2-fold?).

What about the choice between a two- or a one-tailed t-test? Although the one-tailed test is tempting since it is less stringent than the two-tailed test (the threshold for calling a t-statistic significantly different from 0 is lower in a one-sided test), you are only allowed to use a one-tailed test if you know upfront (before you do the experiment) that you are only interested in one-sided differences. For instance, if you make a knock-out, you may assume that the gene will be downregulated in the mutant or not changing expression (if the mutant is not ok) but you do not expect to find overexpression of the gene you knocked out. In this case you could use a one-tailed test to check if the gene is properly knocked out.


Two tails comparison

Select the two-tailed test, with the lowest significance level (0,01). The groups are defined using the mouse.

05-define-groups.png

The comparison is based on the selected test method. The choice of two tails allows to find DE genes, either up- or downregulated.

06-start-test.png

The result is a very long list of DE genes with barplot view on the right confirming the difference in expression between the two groups. However, in this format the list is not really useful

07-inspect-results.png

Luckily, the right menus allow us to do more usefull things with the selected genes

07b-inspect-results-more.png

Handicon.png Filters and Find related data are not detailed here and left to your curiosity

1. Profile data

This lets you download differentially expressed data selected by the tool

07c_download-profile-data.png


07d_download-profile-data_result.png

2. Profile pathways

Search pathways enriched in the list of DE genes

07e_find-pathways.png


07f_FLINK-found-pathways.png

Column legends

07g_column-legend.png


Rank mean difference "4x & either"

The two-tailed t-test found too many DE genes to be specific for given pathways. This is because of the experimental setup, comparing two completely different tissues (heart and diaphragm). In most microarray experiments, the groups that you compare are not that different and you will find less DE genes. The only way we can reduce this number is by doing a test that evaluates fold changes and selecting a high threshold, e.g. the rank mean difference test with a 4-fold difference threshold. Now, we get only the 152 genes with the most extreme differences in expression between heart and diaphragm (at least a 4-fold difference).

Search pathways enriched in the obtained subset

07h_rank-mean-difference.png


07i_found-genes.png

The pathways are now very specific to heart biology as expected

07j_FLINK-found-pathways.png

 


 

Cluster heatmaps

Two options are presented in this part, the hierarchical clustering and the K-Mean approaches.

Hierarchical clustering

Hierarchical clustering and heatmaps are classical representations of differential expression that nicely provide high level view on the data.

You may choose between Euclidean distance and Pearson correlation as measures for the distance between two genes (expression profiles). You must realize that the data are not standardized by this tool so both distance measures will measure different things:

  • Euclidean distance will sum up absolute differences between expression levels: it groups genes that have similar expression VALUES
  • Pearson correlation will measure the similarity of the expression profiles: it groups genes that have similar expression BEHAVIOUR, up in the same samples and down in the same samples

So the two distance measures will result in different trees, as you can see in the figure below showing the immediate neighbours of Creg1 (represented twice on the array) in the resulting tree.


hierarchical.png

You see that Euclidean distance groups genes with very similar expression values (to the same extent high or low), while Pearson correlation focusses on trends (high in the same samples regardless of how high exactly) and that they result in different groupings.

So it's important to make the right choice of distance measure for what you want to know !


08-hierarchical-config.png

You see that you also have to choose how to define the distance between clusters. This choice is not as important as the choice between Euclidean distance and Pearson correlation. The difference between the three proposed methods is:

  • average: you measure the distances between cluster centers, the standard choice
  • complete: you measure the distance between the farthest members of each cluster
  • single: you measure the distance between the closest members of each cluster

A number of interactive links allow changing the correlation method, changing colors, and/or selecting heatmap regions to plot other graphs or export data to file

09-hc-results.png


09b-hc-results-selection.png

 


 

K-Means clustering

.

Technical.png If you run K-means repeatedly, it will always return different clusters. So you will not get the same results as the ones below!!! BUT you will find back mostly the same genes in similar clusters

The user decides arbitrarily the number of clusters to be found. Too few will lead to mixed profiles while too many will lead to apparent redundancy.

10-kmeans-config.png

From the obtained solution, it is clear that the four clusters are grouping genes having clear trends in both sample groups. Groups #1 and #3 being highly contrasted between Diaphragm and Heart while #2 and #4 are not/more subtle DE groups

11-kmeans-results.png

Clicking on the first group to get more details for genes higher in the Diaphragm

12-up-down-cluster.png

Similarly, for genes higher in the Heart

13-down-up-cluster.png

Finally, zooming on a small number of genes allows plotting their profiles and exporting the results to file

14-line-plot.png

Experiment design and value distribution

This QC tool plots box-plot for each sample and allows evaluating the global quality of the dataset. A good dataset has a constant median line and similar distributions.

15-distribution-box-plots.png

download exercise files

Download exercise files here.

Use the right application to open the files present in ex3-files

References:
  1. http://www.ncbi.nlm.nih.gov/sites/GDSbrowser/
  2. http://www.ncbi.nlm.nih.gov/geo/info/datasets.html

[ Main_Page | Hands-on Analysis of public microarray datasets | PubMA_Exercise.2 | PubMA_Exercise.3 |
| PubMA_Exercise.4 ]