A reanalysis of mouse ENCODE comparative gene expression data

Recently, the Mouse ENCODE Consortium reported that comparative gene expression data from human and mouse tend to cluster more by species rather than by tissue. This observation was surprising, as it contradicted much of the comparative gene regulatory data collected previously, as well as the common notion that major developmental pathways are highly conserved across a wide range of species, in particular across mammals. Here we show that the Mouse ENCODE gene expression data were collected using a flawed study design, which confounded sequencing batch (namely, the assignment of samples to sequencing flowcells and lanes) with species. When we account for the batch effect, the corrected comparative gene expression data from human and mouse tend to cluster by tissue, not by species.


Introduction
The mouse ENCODE Consortium has collected multiple types of genomic and functional data in order to better understand the potential utility of the mouse as a model system for biomedical research. To study gene expression levels, the Consortium collected RNA sequencing data from multiple tissues from human and mouse. Their comparative analysis revealed that gene expression patterns tend to support clustering of the data by species, rather than by tissue (Figure 2a in reference 1).
This pattern was confirmed and discussed in greater detail in a companion paper by Lin et al. 2 , which also acknowledged that this observation is somewhat unexpected. Indeed, previous comparative studies reported that gene expression data from human and mouse (and across other species more generally) tend to cluster by tissues, not by species. Lin et al. proposed that previous studies might have been biased in their focus on a few 'specialized' tissues that tend to express the largest number of 'tissue-specific genes', while the overall pattern supports less tissue specificity.
The implications of the observation that human and mouse gene expression data may be clustering by species more than by tissues can be profound. To a large degree, modern biology is built upon the empirical observation that homologous gene regulatory networks establish the identities of homologous cell-types, tissues, and organs across species -the results of Lin et al., if true, challenge these observations and the biological basis of homology. From a more practical perspective, the mouse is arguably the most important animal model for biomedical research. If gene regulation in any mouse tissue is markedly more representative of a general mouse regulatory network than the regulatory network of a corresponding human tissue, this would call into question the utility of the mouse, and perhaps any other non-human animal, as a useful model system for biomedical research.
Here, we present a reanalysis of the mouse ENCODE Consortium comparative RNA sequencing data. We argue that a flaw in their study design raises doubt regarding their conclusions.

RNA-Seq data, genome and gene annotation files
In December 2014 we asked and were kindly provided by the authors of Lin et al. 2 the names of the sequence files used in their comparative analysis. Based on this information we obtained sequence files in FASTQ format (Supplementary Table 1) from the ENCODE project 1 site (https://www.encodeproject.org/; some of the files were only available from early January 2015).
Sequencing study design Based on the sequence identifiers found in the FASTQ files, we reconstructed the sequencing study design used to collect the gene expression data in Lin et al. 2 . The sequence identifier line in a FASTQ file generated from an Illumina sequencing run can take two formats, depending on the version of the Consensus Assessment of Sequence and Variation (CASAVA) pipeline used to generate it. Prior to version 1.8 of this pipeline the sequence identifier line was of the following format (CASAVA v1.7 user guide p.88; downloaded from: http://support.illumina.com/downloads/casava_ software_version_17_user_guide_(15011196_a).html @<machine_id>:<lane>:<tile>:<x_coord>:<y_coord>#<index>/ <read_#> Starting from version 1.8 the sequence identifier line is of the format http://support.illumina.com/help/SequencingAnalysisWorkflow/ Content/Vault/Informatics/Sequencing_Analysis/CASAVA/swSEQ_ mCA_FASTQFiles.htm @<machine_id>:<run number>:<flowcell ID>:<lane>:<tile>:<x-pos>: <y-pos> <read>:<is filtered>:<control number>:<index sequence> Below is a sequence identifier line from the mouse pancreas read1 FASTQ file (sequence identifier lines from the remaining FASTQ files were of similar format): @D4LHBFN1:276:C2HKJACXX:4:1101:3448:12374 1:N:0:AGT-TCC Based on this information we inferred that the FASTQ files were generated by CASAVA version 1.8 or higher. Thus, we could extract from the sequence identifiers the following details that pertain to the sequencing study design: machine identifier, run number, flowcell identifier, and flowcell lane number. We found that the sequencing was performed in five batches, each consisting of a multiplexed single run on a single lane on one of four sequencers (Figure 1; note that two of the batches, composed of human samples only, differed only in their lane number). The design was such that only one batch contained samples from both species. The remaining four batches could be divided into pairs where each of the two batches had a nearly identical tissue composition, but a different species.
Ortholog annotation Following Lin et al. 2 , we used the protein-coding ortholog list generated by the modENCODE and mouse ENCODE consortia 5 . A file containing all orthologs from human, mouse, fly and worm was downloaded from http://compbio.mit.edu/modencode/orthologs/ modencode.common.orth.txt.gz. From this list we extracted 14,744 human-mouse one-to-one ortholog pairs, for which both members were included in the transcript annotation files we used. We note that this number is lower than the ~15,106 ortholog pairs reported in Lin et al. We are not certain of the meaning of the '~' in the report of the number of ortholog pairs analyzed by Lin et al. Nevertheless, we believe that a possible explanation for this disparity is a parsing error. The last two columns of the 'modENCODE ortholog file' represent the number of genes from each species in the ortholog group. One of the steps required to obtain the subset of ortholog groups for analysis is to select those records where the two last columns have a value of 1 (i.e. one-to-one ortholog pairs). We found that if this selection is done through a command line search that does not require that the value in the last column be exactly "1", but rather just begins with "1", then the result is 15,104 putative human-mouse ortholog pairs.
Quality assessment of RNA-Seq data We used the FastQC software v0.10.0 (http://www.bioinformatics. babraham.ac.uk/projects/fastqc/) to assess the quality of the individual FASTQ files (Supplementary Table 2-Supplementary Table 6). We were concerned by evidence for GC content bias and overrepresented sequences. To examine the latter in greater detail, we mapped the sequences overrepresented in at least one sample to the genome of the respective species, using BLAT searches 6 against the hg19 (human) and mm10 (mouse) assemblies at the UCSC genome browser site (http://genome.ucsc.edu/) 6 . We found that in both species many of the overrepresented sequences mapped perfectly to the mitochondrial genome (Supplementary Table 3-Supplementary  Table 6). For the mouse pancreas sample only, we also found many overrepresented sequences mapped to regions with rRNA repeats from the SSU-rRNA_Hsa and LSU-rRNA_Hsa families.
Mapping RNA-Seq reads to genome sequences We mapped the RNA-Seq reads to their respective genomes using Tophat v2.0.11 7 with the following options: "--mate-inner-dist 200" (i.e. inner mate distance is 200nt, based on paired-end reads with length 100nt each and an insert size of 350-450nt ); "--bowtie-n" (i.e. the "-n" option will be used in Bowtie 8 in the initial read mapping stage); "-g 1" (i.e. multi-mapping reads will be excluded from alignment); "-m 1" (i.e. one mismatch is allowed in the anchor region of a spliced alignment); "--library-type fr-firststrand" (the libraries had been constructed using the Illumina TruSeq Stranded mRNA LT Sample Prep Kit 2 ). An exception was the mouse pancreas sample, for which the mapping process stalled consistently at the same stage.
For this sample we used Tophat v1.4.1 8 with the same options as above. Tophat requires a Bowtie 8 index. For human we used the Bowtie index that was packaged with the genome sequence in the file downloaded from the Illumina iGenomes page (http://support. illumina.com/sequencing/sequencing_software/igenome.html). For mouse we built an index using the bowtie-build utility from Bowtie v2.2.1 (v 0.12.7 for the index used with Tophat v1.4.1).

Calculating gene GC content
For each of the two species we used the appropriate GTF file to generate a table, which contains for each gene its ENSEMBL gene identifier its common name, and the GC content of the sequence covered by the union of the gene's transcripts. To this end, we first generated a GTF file where overlapping exons from different transcripts of the same gene were merged into a single "exon" with the same sequence coverage, retaining the association with the gene identifier. Next, we computed the nucleotide content of the exons in this new GTF file using the 'nuc' utility from bedtools v2.17.0 9 . Finally, we computed the GC content for each gene identifier by summing the number of 'G' and 'C' nucleotides in its merged exons and dividing by the sum of counts of unambiguous nucleotides in these exons.

Per-gene FPKM values
We used Cufflinks v2.2.1 10 to compute fragments per kilo base of transcript per million (FPKM) values and aggregate them per gene. The only option used was "--library-type fr-firststrand". For the required transcript annotation file ("-G" parameter) we used the GTF file for the respective species described in the "Genome and gene annotation files" section. We then generated a matrix of 14,744 by 26 FPKM values for each gene (in the ortholog table) and sample. While generating this table we noticed that some of the common gene names were associated with more than one ENSEMBL gene identifier. In some cases we determined that this was due to gene identifiers that have been retired from the ENSEMBL database 3 but were retained in the GTF file (27 and 64 retired identifiers for human and mouse, respectively). These retired identifiers were ignored when constructing the FPKM matrix. For the remaining such cases we incorporated the value from the first appearance of the common name. Per-gene raw fragment counts To compute per gene raw counts from the alignment files produced by Tophat 7 , we used the program featureCounts v1.4.4 11 with the respective species' GTF file specified in the "Genome and gene annotation files" section. For all runs we used the following options: "-p" -indicates that fragments rather than reads should be counted; "-C" -indicates that chimeric fragments will not be included in the summarization process; and "-s 2" -indicates that the pairedends are reversely stranded. We next generated a matrix of 14,744 by 26 raw counts for each gene (in the ortholog table) and sample. Since the output from featureCounts identifies genes by their gene identifier (the ENSEMBL identifier in our case), whereas the ortholog table uses the gene's common name to identify it, we used the GC content table, which contains both these identifiers to match counts to the correct row in the ortholog table. As we did when generating the FPKM matrix, we ignored the values from retired ENSEMBL identifiers, and if there were still multiple identifiers for the same common name, we used the value from the identifier that appeared first.

Results
In this reanalysis effort, we focused solely on the RNA sequencing data that can be mapped to coding regions. Lin et al. 2 reported additional results, related to data on the expression of non-coding transcripts and histone marks. We did not reanalyze these additional data types.
Lin et al. 2 analyzed both previously published and newly collected human and mouse gene expression data. The previously published data consist of RNA sequencing from ENCODE, the Illumina Human BodyMap 2.0, and the Roadmap Epigenomics Mapping Consortium. In these previously collected data sets, human and mouse samples were analyzed by different labs at different times, such that there is a clear batch effect that is confounded with species. Lin et al. 2 clearly explains this limitation of the previously published data. They state that in order to address this issue they focus on the analysis of only the newly collected data -RNA sequencing data of samples from 13 human and mouse tissues that were collected by the same lab, using the same sample processing protocol. We focus our reanalysis study on the same newly collected data set (see Methods).
Replication of sample clustering by species As a first step of our study we set out to replicate the analysis of Lin et al. 2  We log 2 -transformed the FPKM matrix (after adding 1 to avoid undefined values). To visualize the data, we used an approach that is similar in principle to that used by the ENCODE mouse consortium and Lin et al. Specifically, we used the function 'prcomp' (with the 'scale' and 'center' options set to TRUE) to perform principal component analysis (PCA) of the transposed FPKM matrix (so that samples were now in rows and genes in columns), after removal of invariant columns (genes). Scatter plots of the PCA results were generated using the ggplot2 package 13 . In agreement with the findings of Lin et al. 2 the samples cluster mostly by species (Figures 2a, Figure S1 and Figure S2). We also plotted the heatmap of the matrix of Pearson correlations between the 26 samples, using the pheatmap function from the pheatmap package v1.0.2 14 with default settings (i.e. complete linkage hierarchical clustering using the Euclidean distances). Again, samples from the same species tend to cluster together ( Figure 2b). Analysis of normalized data after accounting for batch effects yields clustering by tissue A previous evaluation of normalization methods for RNA-Seq data 15 suggested that FPKM values were not optimal for clustering analysis. Therefore, as a basis for our reanalysis, we used the matrix of per-gene raw fragment counts. The entire analysis was done within R environment v 3.1.3 GUI 1.65 Snow Leopard build (6912) 12 . See Supplementary Text 2 for detailed commands, and a supplement zip file for the R input (available in Zenodo: http:// dx.doi.org/10.5281/zenodo.17606).
Following Li et al. 16 , we removed the 30% of genes with the lowest expression as determined by the sum of fragment counts across all samples. Next, due to the presence of mitochondrial genes among the overrepresented sequences in the data, we also removed reads that map to the 12 mitochondrial genes. This left us with expression data from 10,309 genes for analysis. We note that merely limiting the analysis to this subset of genes does not have a marked effect on the patterns reported by Lin et al. (Figure S3; detailed commands in Supplementary Text 3, and a supplement zip file for the R input (available in Zenodo: http://dx.doi.org/10.5281/zenodo.17606)). We performed within-column normalization to remove the GC bias in the data, indicated by the initial quality assessment. To this end, we applied the 'withinLaneNormalization' function from the EDASeq package v2.0.0 17 to each column in the matrix, using the gene GC values for the species associated with the column. Next, we used the 'calcNormFactors' from the edgeR package v3.8.6 18 , with the trimmed mean of M-values (TMM) method 19 , to calculate normalization factors for the library sizes for the samples. We used these normalization factors in the depth normalization of the columns (using the column sums of the original, unfiltered, counts matrix as a proxy for library sizes). The normalized data were log 2 -transformed (after adding '1' to each value in the matrix to avoid undefined values).
We then considered how to account for the fact that the assignment of samples to sequencing flowcells and lanes was nearly completely confounded with the species annotations of the samples (Figure 1). The consideration of 'batch effect' was the most important difference between the analysis that recapitulated the patterns reported by the mouse ENCODE papers (the previous 'Results' section) and the current reanalysis effort. Specifically, we accounted for the sequencing study design batch effects using the 'ComBat' function from the sva package v3.12.0 20 , with a model that includes effects for batch, species and tissue. For this purpose the samples were classified into five batches, based on the sequencing study design (see methods and Figure 1).
To visualize the data, we used the function 'prcomp' (with the 'scale' and 'center' options set to TRUE) to perform principal component analysis (PCA) of the transposed log-transformed matrix of 'clean' values (after removal of invariant columns, i.e. genes), and the ggplot2 package 13 to generate scatter plots of the PCA results. None of the first five principal components (accounting together for 56% of the variability in the data) support the clustering of the gene expression data by species (Figure 3a and Figure S4- Figure S5). However, the sixth principal component, which accounts for 6% of Figure 3. Clustering of data once batch effects are accounted for. a. Two-dimensional plots of principal components calculated by applying PCA to the transposed matrix of batch-corrected log-transformed normalized fragment counts (from 10,309 orthologous gene pairs that remained after the exclusion steps described in the results) for the 26 samples, after removal of invariant columns (genes). b. Heatmap based on pairwise Pearson correlation of the expression data used in panel a. We used Euclidean distance and complete linkage as distance measure and clustering method, respectively. the variability in the data, does support such a clustering, suggesting that even though the 'species' and 'batch' variables are confounded, accounting for 'batch' does not remove completely the variability due to 'species' ( Figure S5). We also plotted a heatmap of the matrix of Pearson correlations between the 26 samples, using the pheatmap function from the pheatmap package v1.0.2 14 with default settings (i.e. complete linkage hierarchical clustering using the Euclidean distances). This time the heatmap shows considerable clustering of the comparaive gene expression data by tissue ( Figure 3b).

Discussion
In our reanalysis we have made a number of specific choices, including the exclusion of a certain subset of lowly expressed genes, the specific approach we chose to summarize the count data, the standardization and normalization methods we used (for example, we chose to standardize by the total count of reads that mapped to the ortholog gene pairs), the approach we used to account for the GC content bias, and the method we used to account for the sequencing design batch effect. Moreover, we excluded the sequencing data from 12 mitochondrial genes from both species, a step that -to the best of our ability to determine -was not taken by the original studies. In addition, our definition of ortholog gene pairs differs slightly from that of the original study, as we discussed in the methods. In practice, only the correction for the sequencing design batch effect had a drastic impact on the results. For example, without accounting for batch, using per-gene raw fragment counts instead of FPKM values does not seem to impact the degree to which the uncorrected data support clustering by species ( Figure S6).
Visualizing or plotting the data is another important area where different choices can sometime lead to quite distinct conclusions. We chose to display, in addition to the PCA plots, heatmaps based on the correlations among the samples. We note that if the actual data (not pairwise correlations) are clustered, the observed patterns (by species or by tissues, in the respective analyses), seem practically identical ( Figure S7). The heatmaps shown in the main figures are based on Pearson pairwise correlations, which provide the highest level of clustering by tissue in the analysis that takes into account batch effects. Alternative heatmap plots based on either Spearman pairwise correlations or other distance measures and clustering methods look similar in principle ( Figure S8 to Figure S10), but the clustering by tissue is somewhat less pronounced (clustering by species, when batch is not accounted for, is more pronounced).
It is important to note that most of the analysis and plotting decisions we have made contributed to a somewhat better clustering of the expression data by tissue, both visually and empirically. We have made these -mostly standard -analysis and plotting choices regardless of the end result (namely, we believe that these are objectively reasonable choices). Importantly, we made identical choices for the clustering analysis and plot types for the data with and without batch correction, and our conclusions are robust with respect to a wide range of possible alternative approaches ( Figure S7- Figure S10).
That said, we do acknowledge that we find the clustering of the data by tissue to be a more intuitive pattern. In other words, we believe that the clustering of comparative gene expression data by species -a result that contradicts previous observations -is a surprising outcome. Hence, we would have intuitively accepted as more correct most reasonable choices of analysis pipelines and data visualizations that supported a greater degree of clustering by tissue.
As we mentioned above, most of the choices we made resulted in little difference to the overall pattern. It was only the correction for the sequencing design batch effect that had a profound impact. Once we accounted for the batch effect by using ComBat, the comparative gene expression data no longer clustered by species, and instead, we observed a clear tendency for clustering by tissue. This is not surprising, as the sequencing batch, which we corrected for, was nearly entirely confounded with species. It stands to reason that some individual gene expression levels do cluster by species and some by tissue (see for example, Figure S5). While previous data sets strongly support a general clustering of gene regulatory phenotypes by tissue 21 , we expect the degree of clustering of the gene expression data to differ somewhat across tissues. Yet, in this particular case, by removing the confounding sequencing batch effect we also removed most of the species effect on gene expression levels (a similar case of confounding batch and main effect of interest was discussed a few years ago, with respect to gene expression differences between human populations 22 ).
One could potentially employ more sophisticated modeling approaches to try and estimate separately the batch and species effects. One idea would be to rely on the fact that there are five sequencing batches, but only two species. This, however, is complicated by the fact that the two sequence batches specific to the human samples share the same run and flowcell (potentially a smaller batch effect), while the two sequence batches specific to the mouse samples extend over different instruments (potentially a larger batch effect). In any case, we feel that such modeling is beyond the scope of this reanalysis effort. Instead, we conclude that the study design used by the mouse ENCODE consortium was flawed with respect to the questions they set out to address.
In summary, we believe that our reanalysis indicates that the conclusions of the Mouse ENCODE Consortium papers pertaining to the clustering of the comparative gene expression data are unwarranted. In the narrow context of our reanalysis effort, we state that their conclusions are unwarranted, not wrong, because the study design was simply not suitable for addressing the question of 'tissue' vs. 'species' clustering of the gene expression data. That said, a large body of independent previous work supports general clustering of comparative gene expression data by tissue.
Finally, we note that in this reanalysis effort, we have only focused on the RNA sequencing data collected by the mouse ENCODE consortium. We have not considered information with respect to the study design used to collect the many other types of data reported by this consortium. Given our findings, we believe that it is appropriate to call for a careful review of these other data sets as well.

Data availability
All data are available from the Mouse ENCODE consortium; see Table S1 for specific source URLs and accession numbers.

Software availability
We provide supplementary files of the python codes used to process and prepare the data for analysis with R, and the data files for the python codes. We also provide the R codes we used to perform the different analyses as supplementary files, as well as the input for the R codes.

Competing interests
The authors have no conflicts of interest or competing interests to disclose.

Grant information
This work was supported by NIH grant MH084703.
I confirm that the funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Acknowledgment
We thank J. Lieb, V. Lynch, J. Novembre, L. Pachter, D. Graur, M. Stephens, N. Banovich, and J. Leek, for comments on the manuscript. We also thank 5 additional colleagues, who asked us not to reveal their names, for their valuable comments. Table S1. Source of RNA-Seq data. Table S2. Summary of test scores for the 52 FASTQ files analyzed. Table S3. Overrepresented sequences in read1 human files. Table S4. Overrepresented sequences in read2 human files. Table S5. Overrepresented sequences in read1 mouse files.  that expression data from human and mouse cluster more by species than by tissue.

Supplementary tables
The Gilad--Mizrahi-Man paper consists of three "results": A report of the experimental design in Lin et al.

1.
An attempt to reproduce the results of Yue et al. and Lin et al. that pertain to the claim about species vs. tissue clustering of expression data.

2.
A re-analysis of the Lin et al. data in a manner that addresses shortcomings in the original experimental design.

3.
The first is the observation that Lin et al. improperly designed their experiment by confounding species with batches sequenced, thereby leading to a possible "batch effect" affecting their results. This observation was already published as a preprint by the first author on the pre-print server Twitter (see https://twitter.com/Y_Gilad/status/593088451462963202). It is evident that while the Lin et al. results may, in some technical sense, be "reproducible" they were certainly not "usable" as published. Gilad--Mizrahi-Man carefully expose a vast number of choices in software and processing options poorly described in Lin et al., and whose effect on the final result(s) is unclear. To quote just one example, they write that "An exception was the mouse pancreas sample, for which the mapping process stalled consistently at the same stage", a problem that led them to use TopHat v1.4.1 instead of TopHat v2.0.11. One may wonder whether software choices and other decisions in analysis affect final results, and Gilad--Mizrahi-Man address this question (although only partly).
For example, one fundamental analysis choice is whether to quantify abundances of genes by summing raw "fragment counts" from alignments to gene regions, or via the summing of abundances as quantified by probabilistic assignment of ambiguously mapped reads. Gilad--Mizrahi-Man cite a paper by Dillies et al. (and the French StatOmique Consortium) suggesting that "FPKM values were not optimal for clustering analysis" to argue for using "fragment counts". I strongly disagree with this choice because transcript abundances are necessary to accurately estimate gene-level abundances, a point that Dillies et al. fail to realize. As pointed out in my own paper on Cufflinks 2 (Trapnell et al. 2012) wrong does not cancel wrong for differential analysis, nor does it for the purpose of clustering.
In any case, Gilad--Mizrahi-Man do examine whether quantification by EM affects results and in a later statement they state that "using per-gene raw fragment counts instead of FPKM values does not seem to impact the degree to which the uncorrected data support clustering by species", a result summarized in their Figure S6. While I applaud them for checking the dependence of results on this choice, without further analysis the question remains of whether other analysis choices affect results (although to be fair to the authors, the number of tests that would have to be conducted is enormous and quite possibly practically intractable). Nevertheless, it would be interesting if, for the purpose of future transcriptomics analyses, Gilad--Mizrahi-Man were to investigate some key steps as to their effect (e.g. annotation, an issue discussed recently by Ongen and Dermitzakis, or mapper choice).
The final result of Gilad--Mizrahi-Man is a re-analysis of the Lin et al. data from which they observe that a basic correction for batch effect removes the strong clustering of expression profiles by species touted in Yue et al. and Lin et al. The question of the relative species/tissue contribution to expression profile is of course fundamental and interesting, and obviously further data, carefully curated and analyzed, will answer the question definitively. As far as the Lin et al. paper goes, the Gilad--Mizrahi-Man paper certainly casts doubt on the suitability of the data for answering the question. For one thing, the term "batch effect" is unfortunately rather generic and in this specific case that has become a problem. After initial posting of the preprint by Gilad on Twitter, the authors of Lin et al. resequenced their libraries in a different configuration, but further investigation of the experimental design by Gilad et al. (subsequent to initial posting of a preprint of the article I'm reviewing) appears to have revealed additional problems. For example, tissues in human and mouse were selected from males vs. females respectively (except in ovaries and testis) resulting in another potential bias that would skew expression profiles to cluster by species rather than tissue. In other words, with their paper and subsequent analysis Gilad and Mizrahi-Man have convinced me to be skeptical of the data and conclusions of Lin et al. But that is not really the point. What matters now are the carefully documented serious shortcomings in the computational and experimental methodology of Lin et al. and for this reason I have approved the Gilad--Mizrahi-Man manuscript. Hopefully the issues raised will be properly

Mick Watson
The Roslin Institute, University of Edinburgh, Edinburgh, UK The study is carried out well and the results support the conclusion.
The paper would benefit from including some of the discussion points made in the comments made to v1 of the paper. Lin et al. have re-sequenced the samples and removed the sequence lane batch effect, and reproduced the same result; however, the samples themselves are confounded, in that they were treated differently prior to sequencing. This discussion should be added to the paper.
I would also like to see a discussion of artifacts which are discoverable within the data, for example: the human samples have significant numbers of rRNA reads compared to the mouse samples. This should not happen with mRNA-Seq which includes a polyA selection.

expression-units/
Competing Interests: It may be apparent from the review, but we have also been re-analysing these data.
In this paper Gilad and Mizrahi-Man reanalyze a high-profile dataset from the ENCODE consortium that was used to argue that gene expression levels are more similar for different tissues from the same species than the same tissues from different species, a somewhat counterintuitive result that contradicts earlier claims (including those by Gilad).
The main result of this new work is a simple observation: in the original experiment described in Lin et al., samples from the same species were run in the same sequencing batch. Since there are well-known batch effects, this is poor experimental design that calls into question the claim by Lin et al. that data clusters by species, since the data could instead by clustering by sequencing batch.
There is always a bit of a challenge in figuring out what to do with an observation. As the authors point out, this aspect of the experimental design effectively renders the data useless for asking questions about the relative contribution of species and tissue to gene expression variation. Yet it would seem like too light of a paper to simply say "The original authors messed up. Their claims are therefore invalid. QED." So this paper contains a few analyses designed to ILLUSTRATE the point they make. They are not exactly results, since, once you realize that batch and species are completely confounded, correcting for batch will inevitably remove the species signal. In this context, it's a bit weird to present such an analysis as if it is a result, but I don't really see a way around it, and the authors are forthright in pointing out that their main observation is that the data from Lin et al are useless for addressing this issue, and that they can make no claims, even after correcting for batch effects, about what the data actually do say.

Competing Interests:
No competing interests were disclosed.
I confirm that I have read this submission and believe that I have an appropriate level of expertise to confirm that it is of an acceptable scientific standard.

1.
As mentioned, with the existing study design it is impossible to completely tease out species from batch. However, there is a relatively simple data analysis that can be performed to explore the possibility that instrument or lane are a large enough source of variability to overcome the tissue effect, which is know to be large. The analysis is: i) perform the same PCA analysis on the mouse data and compare the two instruments and then 2.
ii) perform PCA analysis on the human data and compare the two lanes.
If in fact lane and instrument are a large sources of variability we should see it here. Of course, there is still the possibility that the instruments used for humans was very different to the one used for mice, while the two instruments used for mice were similar. Due to confounding we won't know for sure, but the analysis described here will at least give us at least a lower bound on how large these effects can be.
There is a comment in the F1000Research article from the first author of the PNAS article describing a second experiment in which confounding with instrument or lane was not present. In this analysis species continues to be the first few PCs. In a second version of this article, the authors can perhaps comment on this, as well as some of the other comments that suggest other possible sources of variability that may be confounded with species. Thank you for spending the time to provide a review of our work. We agree with you that given the study design used by the mouse ENCODE consortium, applying a batch correction is futile. Indeed, we explicitly explain that in our discussion (you referred to that section of the text in your review).
We further agree that it would be intellectually interesting to research the extent of the batch effect further -for example, by following your suggestion on how to test for the effect of instrument and lane.
However, we feel that this additional effort is beyond the scope of our study. The mouse ENCODE consortium papers did not discuss (or account for) the sequencing study design. We spent considerable effort tracking the details that allowed us to reconstruct their design. We pointed out in our paper that given this study design, the unusual biological result reported by the mouse ENCODE consortium might have a technical explanation. We believe it is the responsibility of the mouse ENCODE consortium authors to provide evidence that excludes this technical possibility, rather than us having to prove that it is indeed the likely explanation.
Which leads us to your third point: Indeed, the mouse ENCODE consortium authors commented that they have now collected additional sequence data, using a different design, and that their results held. In that sense, we believe that this means that the mouse ENCODE consortium authors accepted our claim that their original design was flawed.
Yet, as mentioned in a few other comments here, there is an additional technical batch effect that was not yet excluded -related to tissue extraction and sample preparation. We plan to discuss this additional technical batch effect in a revised version of the text (we will wait to see additional reviews before we provide a revised version of the paper).
Again, thank you for your time and thoughts.
Competing Interests: No competing interests were disclosed.

Comments on this article Version 1
Reader Comment 07 Aug 2015 Lenny Teytelman, Protocols.io, USA There continues to be a lot of discussion about the proper way to review and comment on research post-publication. Many believe that raising questions in public without contacting the journal or authors first should be a prerequisite. I disagree.
There are many cases where contacting the original authors for clarification and discussion makes sense and should be done. However, voicing concerns and disputing conclusions and results that have been published can not be conditional on the approval of the author, editor, journal. This has never been the case in science.
We Our principal claim is NOT that the data cluster by tissue rather than by species. Rather, our claim is that your study design does not allow you to address the question of whether the data cluster better by tissue or species. Yes, it is true that my own intuition is that your conclusions are wrong, but using your data I cannot provide strong support for my intuition -and you cannot support yours. Your study design is simply flawed with respect to this particular question.
Let me address directly a few of your points to further explain my position: Clustering analysis (your third point): Indeed, it is true that choosing different approaches to cluster the data will result in different patterns. But that is exactly my point! Your observations are not robust. The choices Orna and I made are not unusual. Moreover, as we have shown in the large ppt we provided ("the 80 slide opus" as Mike referred to it) -we have seen similar results using a range of standard approaches for clustering. We do not claim that it is impossible to achieve better clustering of the data by species by using a specific approach -rather, we claim that since the results depend on the particular choice of clustering approach, they are not robust and therefore the conclusions are not well supported.
Linear modeling (your fourth point): I have had this discussion with Mike at Banff, and he conceded that this is indeed a problematic issue in your analysis (please see the recording of the Q&A where we had this exchange at http://www.birs.ca/events/2015/5-dayworkshops/15w5142/videos/watch/201508040909-Snyder.html). Briefly, the problem is that you only have a single observation from each tissue-species combination. Beyond the fact that I am not sure how to implement an effective optimization of a mixed model with three random effects, your error, which you also need to estimate, is effectively confounded with the interaction effect of species by tissue. Yet, because there are no replicates -you simply do not have enough degrees of freedom to estimate that. I think (but I am not sure) that this is the reason that the residual is so much larger than any other effect. In other words, using your own model, you report that most of the data is simply noise… That is actually not true (I hope you agree with that…) and is an artifact of using the wrong model given the study design.
Lastly with respect to sex and age effects (your fifth point): This part of the comment demonstrates well my claim regarding the flawed study design. Your design does not allow you to directly model and estimate all of these effects in the context of tissue and species effects because everything is confounded. Because of this, you have to use your and/or other data to try and independently estimate the effect size of the different covariates and compare them to estimates of the effect sizes of the confounded tissue and species effects. But this is a flawed way to analyze the data. This is exactly why I claim that your design is flawed. You should have been able to directly account for these effects in your own data, by using an effective study design. As is, we can continue to argue about our intuitions regarding tissue and species effects, but the data you collected is simply is not appropriate to actually test these notions. II. PC1 showed tissue-specific clustering; PC2, species-specific clustering. III. With the re-sequenced data, Gilad and colleagues showed clustering which was mostly tissuespecific. IV. Tissue effects accounted for more variance than that of species. V. The mouse and human samples of Lin et al. were not matched for age and sex.

I. Turning of the axes
To clarify the relative distribution of variance in the principal components, we have made 2D plots ( https://www.dropbox.com/s/7r7huteeg3qly5o/Fig1.tif?dl=0). In both the original and re-sequenced data, species-specific clustering is apparent in PC2 and 3; tissue-specific clustering, in the lower order principal components. Thus, the general patterns of the original and re-sequenced data have not changed.
In contrast, hierarchical clustering on the first 3 PCs shows strong species-specific clustering (right plot https://www.dropbox.com/s/14i42swag9fwqcx/Fig3.tif?dl=0). One of the main observations of the Lin et al. paper is that samples in the top principal components cluster largely by species. This finding holds for both the original and re-sequenced data as well as the combined data from five different labs (Fig. 1A of Lin et al. 1 ). Additionally, we have replicated the pattern on human and mouse CAGE-seq data from the Fantom5 project 2 (see previous F1000Research and Pubmed posts).

III. Clustering of the re-sequencing data
The RNA-seq data of the original and re-sequenced sets should not be expected to be identical for several reasons: The second set did not include ovary (exhausted material), different model sequencing instruments were used, and the reagents from the supplier were different. Also, the RNA may have changed over the two years since their extraction from tissue, especially for spleen and pancreas. Finally, the multiplexing design of the re-sequenced data favors tissue clustering.
With three of the thirteen tissues eliminated from analysis of data under the new multiplexing design, Gilad et al. were able to perform hierarchical clustering with seven tissues matching. We note that there are many ways to cluster data. Depending on the specific method chosen for the more limited data set, clustering can be shown to be mainly by tissue or species ( https://www.dropbox.com/s/ah9wlpigg0wevu7/Fig4.tif?dl=0); yet others show a blend of tissue and species-specific clustering.
Our contention continues to be that tissue composition heavily influences the final clustering pattern. We have shown, as with other investigators 3,4 , that using only five tissues with many tissue-specific genes leads to tissue clustering (Fig 1F of Lin et al.). When we include many more tissues, the clustering is by species (Fig 1A and C in Lin et al.). When the number of tissues is intermediate, the clustering is blended or may go one way or the other depending on method of clustering.  .12), and, species (1 and 2), respectively. y ijk is the log expression value; g i , the gene random effect; t ij , the tissue random effect; s ik , the species random effect. Though the tissue and species effects are nested within that of gene, the relationship between the former two effects are in the context of a crossed design, so additional replicates are not necessary.

IV. Species versus tissue variance
The rationale for the random effects model is as follows. Take the species effect s ij , which must, of course, account for human and mouse and thus require subscript k. However, it must be different for each gene, because for some genes, the difference will be positive, others negative, and the remainder near zero (thus requiring subscript i). The estimate of the actual species effect for each gene is not of interest; instead, the focus is on the variance of the distribution from which it is drawn. Thus, these terms are best modeled as a random effect rather than fixed effects. The same reasoning may be applied to the tissue effect.
We have fit model 1 using library nlme 5 in R (see Appendix for code https://www.dropbox.com/s/axb91137rdyw6x2/Appendix.pdf?dl=0). Analogous to what was seen in the PCA plots, analysis only of tissues expressing high numbers of tissues-specific genes (testes, brain, heart, liver, and kidney) increases the tissue variance component. This finding supports our claim that tissue composition alters the relative effect of species and tissue.
Of note, variance components are the square error divided by degrees of freedom. Therefore, for the five tissue analysis, tissue accounts for more variance than species overall, but after dividing by degrees of freedom, the species variance component is still greater than that of tissue.

V. Effect of sex and age
To examine the impact of sex and age, we performed PCA of GTEx tissues (these data had already been analyzed for different purposes in the paper Lin et al. 1 ) and found that the predominant clustering pattern distinguishes between neural and non-neural tissues ( https://www.dropbox.com/s/3n1tpfa94zkjyia/Fig5.tif?dl=0). Neither sex nor age contributed to the clustering in the first two principal components. Thus, these effects are likely to be small and have not been included in the studies of others as well.

V. Conclusions
We have performed new experiments and analyses to address questions that have been raised about our work by Gilad et al. as well as other researchers in the field. These results as well as analyses of others' data (i.e. GTEx and Fantom5) appear to support our earlier conclusions. These observations hold for a variety of analytical methods-PCA, hierarchical clustering, and regression modeling. Use of these methods to compare the relative effects of species and tissues on transcription have been used in a similar context in other studies 4 . As with others, we did not account for other potential confounders such as GC content differences between species.
In the case in which our results differ the most with those of Gilad et al.'s (viz. clustering of the resequenced Stanford data set minus spleen and pancreas), we explain their increased tissue clustering to be due to their particular clustering parameters as well as the intermediate numbers of tissues examined. Even in this case, however, the partition of variance by the random effects model reveals the species variance component to be greater than that of tissue.
In both the original and re-sequenced data, we observe tissue-specific clustering using data from various sources when only examining a small subset of these tissues-testes, brain, heart, liver, and kidney (Fig. 1F of Lin et al. 1 ). The clustering becomes more species-specific when we include a broader set of tissues. The composition of tissues in the analyses determines the relative effects of species and tissues; consideration of more tissues brings the species differences to the fore. This has been our primary contention for both the original and re-sequenced data analyses.
These observations indicate to us that significant differences exist between the transcriptomes of mouse and human tissues at the level of transcriptional regulation. These relationships are not simple, however, and vary by tissue. Beyond the first-step analysis of clustering relationships, we explore these issues in greater detail in Lin et al. Now more to the point: We agree that the results from PCA are not very insightful in this case. We provide PCA plots because it is the framework of analysis used by the mouse ENCODE authors. Our goal was to show that even using this framework, we couldn't recapitulate the (somewhat incoherent) basis for their biological conclusions. We could try to perform the analysis you suggested based on clustering of samples using the top N PCs, but would the results truly be more insightful? As you pointed out -and we agree -PCA is not really the right tool to answer these questions.
So we might turn to ANOVA. As you wrote, this is a much more intuitive tool with which to address the question. Yet, analyzing all genes together makes little sense, we -again -agree. It does not take into account differences between genes. In our mind this 'combined analysis' actually makes as much sense as using PCA… It provides some vague estimate of the proportion of variance explained by 'species' and 'tissue' for the entire data set. For whatever its worth, we find that species explains very little in this type of analysis (in the new data), while 'tissue' and 'tissue by species' explain a bit more. Note that in the slides we reported the % of 'explained variance' (in fact, most of the variance, in either the old or new data, remains unexplained). This may have been unclear or even a bit misleading and we will change it. Yet -again -there is certainly no evidence for robust clustering by 'species' (Actually, in either the 'old' or 'new' data).
What we really need, is a gene-based analysis. However, the study design makes it difficult to effectively use ANOVA in this case. The analysis by gene is hopeless. With exactly one sample from each variety (one sample from each tissue in each species), one can't estimate a gene-specific interaction effect, one can't effectively estimate the error term, and tissue and species are not always orthogonal either. Combine this with the fact that tissue samples came from different (and unbalanced) individuals, sexes, and ages, and it's a nightmare for analysis. For whatever its worth, we have performed this analysis anyway, ignoring all of the obvious caveats. We will add the slides to the ppt. You will see that, in the new data, a higher proportion of variation in many more genes is explained by 'tissue'.
To us, the most important aspect of your comment is that you raise additional questions about the analysis (both ours and the original mouse ENCODE paper). We believe that everyone will intuitively agree that 'some' genes are expressed more similarly across species and 'some' genes are expressed more similarly across tissues. The conclusion of the mouse ENCODE paper, however, was that "…in general, differences dominate similarities between the two species." (quote from the abstract of the mouse ENCODE PNAS paper).
Yet, the study design, as we stated multiple times, is flawed. One simply can't effectively address the question of the relative contribution of 'tissue' vs' 'species' using this flawed design. That is our only conclusions.
Do you not agree, based on our observations, this discussion, as well as your own concerns about the framework of the analysis, that the conclusion of the mouse ENCODE paper is unwarranted? There are 2 issues that I believe this analysis tries to address (Qi) Is the new data different from the older data.
(Qii) Does tissue or species dominate the variation in new data Regarding (Qi), I think there is no debate that the data are (and should be) different. I believe the PNAS paper authors also agree with this statement. (Qii) is where I think there is still disagreement between the PNAS authors and Gilad et al.
We have atleast 3 main factors contributing to gene expression variation in this dataset -genes, species and tissues. And then several other confounding factors e.g. age, sex. Lets ignore these for now since the analyses in the slides don't directly model these.
(1) The PCA is decomposing the correlation/normalized covariance across the samples featurized by genes via an orthogonal transformation. What we obtain is a projection of the samples into a new orthogonal subspace of PCs (metagenes). Whats nice about PCA is that a small number of PCs (5-10) can potentially explain a large proportion(>80-90%) the covariance captured by 1000s of features (genes) in the original space. This dimensionality reduction allows (arguably) interpretable visualization of the samples relative to each other and can also potentially (not always) help get around the curse of dimensionality that plagues analysis of high-dimensional data (more on this in the next few paragraphs).
I very much like the visualization of the projection of the samples on to each PC as it gives a nice intuitive feel for what the PCs are doing. E.g. Slide 10 shows projections of the new samples on PC1. Several matched tissues from human and mouse project onto the same point (or very close to each other) on PC1 with a few notable exceptions i.e. testis, spleena and pancreas. So one could surmise that PC1 is representing some metagene signature that is largely tissue-associated but does have a species-related component as well (due to the exceptions). Slide 11 shows projections on PC2. Wrt. PC2 the samples from the same species can be seen to be closer to each other and one can obtain a clean separation of species. So PC2 can be interpreted to represent some metagene signature that is species-associated. This may lead one to conclude that since PC2 and not PC1 (which explains more variability than PC2) is more associated with separating species then maybe species < tissue. However, IMHO I don't think one can conclude that especially since PC1 and PC2 only capture 20% and 13% of the variability. So there is still a huge amount of variability to be explained. PC3 for example is also more species associated. PC4 is more tissue-associated and so on. This also makes it rather dangerous to visualize the samples using only the 2 or 3 PCs at a time as any conclusion from such a visualization is incomplete. If we really want to figure out whether samples from the same species or tissue are closer to each other (cluster together) in the subspace of PCs, we should be using the top N PCs that explain a large proportion (90-95%) of variability. Then compute explicit distances (euclidean distance would be very justifiable in this orthogonal subspace) between samples in this reduced subspace. One could then compute within species/tissue vs. between species/tissue distances or cluster the samples using the distance matrix computed in this space of PCs. Neither the PNAS paper nor Gilad et al. do this.
In slide 19 and 20, Gilad et al. present two types of clustering. In slide 19, they first compute Pearson correlation between all pairs of samples using all genes. Then they use the correlation values of each sample with all other samples as new features to compute a euclidean distance between pairs of samples. They use complete linkage hierarchical clustering with this distance matrix to cluster the sames. This clustering shows the samples largely cluster by species with a few exceptions. In slide 20, they compute a distance matrix (correlation?) across samples using all genes and then use this directly as a distance matrix in complete linkage clustering. In this version of clustering, there is a sorta 50-50 split between species and tissues. IMHO, neither of these clustering approaches have a reasonable justification. There is significant literature explaining the curse of dimensionality in high-dimensional clustering and a distance measure based on 1000s of genes (with lots of correlated structure between them) is almost meaningless. This wikipedia article summarizes some of the issues https://en.wikipedia.org/wiki/Clustering_high-dimensional_data (Sorry I didn't have time to find a real reference!). Moreover, Pearson correlation is not really a distance metric (although it is routinely used in clustering) and more importantly it focuses very strongly on highly expressed genes. I think Rafa Irizzary has also highlighted this issue in his review. Lior Pachter has written about this issue in a blog post as well ( https://liorpachter.wordpress.com/tag/mahalanobis-distance/ ) with several suggested alternatives such as the mahalanobis distance. As I mentioned in the previous paragraph, since we are in PCA land, why not cluster the samples projected onto the small set of PCs that explain most of the variability. The mahalanobis distance goes one set further adjusting for the variance explained by each PC. Either of these approaches would potentially get around the curse of dimensionality in that atleast the distances are meaningful. There are of course caveats here as as well. E.g. PCA projects the samples in a single subspace and clusters could exist in different subspaces. But it would certainly be better than clustering based on distances computed using all genes as features and directly link the clustering with the PCA.
Finally, personally I don't like the use of PCA (in the PNAS paper or this analysis) to answer the key question here -what is the relative contribution of species and tissue to variance of expression across all genes. It is incapable of explicitly answering this question.
Which brings me to the ANOVA analysis.
(2) The ANOVA analysis is attempting to directly partition the variance in the expression data as a function of species and tissue. However, the figures in the slides don't make much sense to me.
Slide 18 would lead us to conclude that tissue accounts for 94.8% of the total variance. This to me is impossible. I am assuming here that the linear model being used is expression ~ species + tissue + species:tissue. Yoav hinted on twitter that the reason this looks odd is that the interaction term accounts for a significant proportion of variance. If that is the case, it indicates a dependence between species and tissue and one cannot simply look at the main effects to conclude species > tissue or not. Species may be heavily contributing through the interaction term.
However, I have a bigger issue with the ANOVA analysis. If I understood it correctly, this ANOVA model assumes a single linear model for expression variation across all 3 factors i.e. genes, species and tissues .. but uses only 2 of the factors in the model. It makes little sense to me to have a single linear model for all genes across all tissues and species with no term for gene identity. There is no way to explain differences between genes. It is also assuming fixed effects (I assume) which I don't believe is appropriate for this analysis. IMHO, a random effects model with an explicit term for gene identity is far more appropriate. An alternative is the variance partitioning analysis using random effects for species and tissues shown in Fig 2B. of the mouseENCODE Nature paper for each gene separately. If the analysis is on each gene separately then only having a species and tissue term is reasonable. But across all genes, I don't see how this ANOVA model makes sense. I may be misunderstanding something.
I look forward to your comments.
Competing Interests: I am part of the ENCODE consortium and a co-author on the mouseENCODE Nature paper (but not on the PNAS paper). Prof. Snyder and I are colleagues in the Genetics Department at Stanford University and active collaborators. However, these views are entirely my own and do not represent the views of any other members of consortia, organizations or papers that I am part of.

Author Response 25 Jun 2015
Yoav Gilad, University of Chicago, Chicago, USA For those with little time to read the entire comment, and those who are not invested enough in this area to study all the details of the papers and back and forth discussions, here is a quick summary: Mouse ENCODE data were collected using a flawed sequencing design. We uncovered the flaws and discussed the issue in our paper.
• In response, Mouse ENCODE collected new data from the same samples using a corrected sequencing design and reported that their conclusions are unchanged.
• However, we have now found a clear difference between the original and new data, which was not mentioned by the Mouse ENCODE authors.
• Now in more detail: Shin et al collected new sequencing data from the same samples they used in the original study. In contrast to the original gene expression data, the new data were collected using a sequencing study design that does not confound species with sample assignment to lanes and flowcells (see table included in the comment by Shin). The mouse ENCODE authors wrote that after analyzing the new data, they "arrive at species-specific clustering as previously reported." Yet, the figure provided by Shin is a 3D PCA plot, rotated in a way that makes it difficult to see PC1.
We have performed an analysis of the new data as well. Our findings are not consistent with the statement made by the mouse ENCODE authors.
Specifically, we found that in the new data, gene expression levels in tissues from human and mouse cluster significantly better than in the old data. There are a number of ways to show this, but the conclusion is robust. For example, if one chooses the approach used by the mouse ENCODE authors (using fpkm values and examining the clustering of samples using PCA plots), we can visualize the difference between the original and new sequencing data by considering the correlation between the PC1 values of the corresponding mouse and human tissue samples. I will post the figures on twitter (@Y_Gilad), but here are the numbers: For the original data (where batch effect is confounded with species) the Pearson correlation between the PC1 values of the corresponding human and mouse tissues is 0.40 (P = 0.18).
For the new data (collected using a new sequencing design), the Pearson correlation between the PC1 values of the corresponding human and mouse tissues is 0.64 (P = 0.02). Moreover, two samples (pancreas and spleen) have a number of severe and obvious QC issues (for example, the mapping % is much lower than in all other samples). If one excludes those samples, the Pearson correlation between the PC1 values of the corresponding human and mouse tissues in the new data is 0.89 (P = 0.0006).
As I mentioned, this finding -of differences in clustering properties between the original and new data -are robust with respect to choices in methodology and plotting. In other words, anyone who wishes to download and analyze the original and new data can easily observe the differences between the data sets (for example, whether one uses fpkm or raw counts, normalize the data in different ways or not at all, correct for GC content or not, and whether one plots PCA results or heatmaps, etc…).
The difference between the data sets is due to the effect of sequencing batch on the patterns observed with the original gene expression data. In the new data, species clustering is weaker and tissue clustering is stronger.
I have discussed the difference between the two data sets with Mike Snyder and all other coauthors of the Shin et al. paper before I posted this comment. I did this because it was strange to me that Shin did not comment on the difference between the data sets. Indeed, based on Shin's comment one might assume that Orna and I were wrong, and sequence batch does not actually have a noticeable effect on the data. Yet, as I have shown, the effect of sequence batch is evident by comparing clustering patterns between the original and new data.
By discussing this issue with the mouse ENCODE authors I have learned that they acknowledge that a difference between the two data sets exist, but they did not feel that it warranted an explicit discussion. The authors also declined to post an additional comment clarifying this. They wrote to me: "the point of collecting the new data was to see if the clustering conclusion still stands when samples from different species were mixed in the same lane and in our opinion it did." We agreed to disagree. You must therefore see why it is reasonable to hold you to a high standard and require that no other (e.g., technical) considerations can provide alternative explanation to your observed patterns.
In the case of the sequencing study design, you have confounded species with sample assignments. You might have assumed that this type of batch effect is minor, but many othersincluding GTEx -have shown that it is actually quite considerable.
Thus, regardless of your assumptions (or intuition), the confounding sequencing batch effects could have potentially explained your original observations. I think that you have recognized that there is no way to argue against this rationale (based on the original data), which is why you sequenced the samples again, using a different study design.
The problem is that through the Twitter discussion, our colleagues raised an alternative, possibly even more significant, technical explanation for your observations. Based on the details you provide in your methods, it seems that the human and mouse samples were collected using quite distinct protocols. In addition, based on my own experience (and given the description of the tissue collection protocols), I am guessing that the quality of RNA extracted from the human tissues is significantly different from the RNA quality of the mouse tissues (you have not provided RIN data).
So, once again, we suspect a technical confounding batch effect that could potentially explain your observations. One cannot simply assume that batch effects do not exist, or that they are minor (without explicitly testing for them).
Indeed, in our world as scientists, technical explanations must always be excluded or they remain reasonable alternatives. This is especially the case when one's observations contradict those of previous studies. In such cases, in particular, we have a strong prior to favor technical explanations for the discrepancy in observations. It is the authors' responsibility to exclude all of those possible alternative explanations by providing the relevant data.
The recognition of global differences between the human and mouse transcriptomes is consistent with the experiences that many experimentalists have using the mouse model. Given the substantial differences in size and metabolic rates, we do not feel it is implausible that there are strong expression differences reflected in basic metabolic and cellular processes at the organism level. Rather than questioning the utility of the mouse model, which will assuredly continue as an invaluable tool, we propose that a better understanding of these differences between human and mouse may allow us to better utilize this model system as it pertains to the investigation of human diseases.
(***Data mentioned herein will be available for download at the Mouse ENCODE website shortly.) If this is the case, then you do have to normalize by GC content and also gene length before clustering as different degrees of degradation in the samples would lead to consistent loss of RNA from transcripts with specific characteristics across all tissues in the more degraded group.
underwent extensive discussions over a two-year period, and conducted analyses and generated additional data to address concerns pertaining to potential laboratory and batch effects. In an examination of a subset of our data, Mizrahi-Man et al. revisit observations previously encountered by the consortium.
Why then, does the normalization performed by Mizrahi-Man et al. result in tissue-specific clustering? If one normalizes away the species-specific differences, of course one will not see them. In the course of its analyses, the consortium demonstrated that if the mouse and human datasets were separately normalized, the global expression comparisons resulted in tissue-specific clustering (see Extended Data Fig. 1c of Yue et al 2 ). This made intuitive sense to us, as it effectively removed the overall expression differences between the species and made apparent the expression differences between tissues. Indeed, the normalization procedure performed by Mizrahi-Man et al. did just that. Because the sequencing libraries were multiplexed largely by species, their normalization was equivalent to intra-species normalization, which effectively removed the global differences between human and mouse gene expression. To reiterate, this normalization sequence and resultant patterns of data clustering were well known to us and detailed in the ENCODE consortium main paper (see Extended Data There remains the issue of our study design with respect to confounding of lane effect and species. It should be noted that our study design minimized library preparation and primer index effect. A recent GEUVADIS consortium study showed that both factors are each contributors to RNA-seq variance and of much greater effect than that of lane (see Fig. 3c