StatQuest: DESeq2, part 1, Library Normalization

Key Takeaways

stack Quest stack Quest oo I love you stack Quest hello and welcome to stat Quest stat Quest is brought to you by the friendly folks in the genetics department at the University of North Carolina at Chapel Hill today we're going to talk about de seek 2 which is a program people use to identify differential gene expression de seek 2 is a big and complicated program so we're going to break it down into Parts this is part one Library normalization remember rpkm fpkm and TPM those nice methods for adjusting for differences and overall read counts among libraries de 2 doesn't use those methods and neither does EDR by the way why not there are two main problems in library normalization so let's talk about them problem number one adjusting for differences in library sizes for the sake of keeping the example simple let's assume there are only six genes in the genome here we have sample number one which has a total of 635 re reads mapped to it and here we see how the reads are distributed among the six genes sample number two has 1,270 reads mapping to it and we see how those reads are distributed among the genes in the genome sample number one has half as many reads as sample number two the read counts for each gene in Sample number two are twice the read counts in Sample number one this difference is not due to biology but to sequencing depth rpkm fpkm TPM and CPM all deal with this no big deal however there is another problem problem number two adjusting for differences in library composition RNA seek and other high throughput sequencing is often used to compare one tissue type to another for example liver versus spleen it could be that there are a lot of liver specific genes transcribed in liver but not in the spleen this is an example of a difference in library composition you can also imagine seeing differences in library composition in the same tissue type if you knock out a transcription factor or something that regulates Gene transcription let's look at a specific example in this example both libraries are the same size both have 635 reads now assume expression of all genes is the same with one exception assume that only sample number one transcribes the gene called a2m to make matters worse sample number one transcribes a2m at a very high level here we see that 563 reads of a total of 635 reads mapped to this Gene in Sample number one this means that the 563 reads used up by a2m in Sample number one will be distributed to other genes in Sample number two here we see the read counts for all of the genes in Sample number one and Sample number two the read counts for everything but a2m are crazy high in Sample number two however the only differentially expressed Gene is a2m because sample number two does not transcribe a2m all of the other genes get the read counts that would have gone to it and this makes those read counts larger the folks that wrote Deek 2 and EDR were aware that their tools would be used with all kinds of data sets so they wanted their normalization to to handle one differences in library sizes and two differences in library composition we'll start with a small data set to illustrate how DEC 2 scales the different samples the goal is to calculate a scaling factor for each sample the scaling Factor has to take read depth and Library composition into account so the first thing that de 2 does is it takes the log of all the values de seek 2 uses the log base e so these numbers are what we would need to raise E2 in order to get the original value so if the original read count was 10 we'd have to raise e to 2.3 to get that value de seek 2 could have used log base 2 or log base 10 but log based e is the default and R which is the programming language that was used to create de2 anyways I think because it's the default that's why they chose it notice that the log of0 equals negative Infinity this is just because R defines log of zero to be negative Infinity if you'd like to learn more about logs check out the stat Quest on logs the next thing that de seek 2 does is it averages each row anytime you add a number to infinity or negative Infinity you end up with infinity or negative Infinity which is why the average for Gene one is negative Infinity one cool thing about the average of log values is that the average is not easily swayed by outliers to see this let's calculate the average read count for Gene 3 we see that the read counts for Gene 3 in Sample number three are really high that makes it an outlier if we just average the raw read counts for Gene 3 we get 96 now convert the average log value for gene3 into a normal number remember that logs are exponents and in this case they are exponents of e so we have to raise e by 4.3 to get a normal number e raised by 4.3 equal 73.7 the average calculated with the logs is smaller and thus not swayed as much by the outlier and for all you stat questers out there that can remember the names of things averages calculated with logs are called geometric averages hooray we've made it all the way through step two only five more steps to go step three is is an easy step filter out genes with infinity as their average so in this case we're going to filter out Gene number one in general this step filters out genes with zero read counts in one or more samples if you are comparing liver and spleen this will remove all of the genes only transcribed in liver or spleen in theory this helps Focus the scaling factors on the housekeeping Gene genes transcribed at similar levels regardless of tissue type step four subtract the average log value from the log of the counts here we have the log of the counts for each gene in each sample and here we have the average of the log values and all we have to do is subtract that average from each sample so in this case for Gene number two in example number one we subtract 1.7 from 0.7 that gives us -1 for Gene number three in Sample number one we subtract 4.3 from 3.5 that gives us 0.8 and we just do the same thing for all the other genes remember when we subtract the log of one value from the log of another value that's the same thing as the log of dividing those two values or the ratio of those two values so we're really checking out the ratio of the reads in each sample to the average across all samples this will allow us to identify genes within each sample that are expressed at levels significantly higher than the average or close to the average or significantly less than the average okay now let's move on to step five calculate the median of the ratios for each sample so here we have the log of the ratios of the reads for each gene divided by the average for each gene and all we have to do is calculate the median value for each sample note using the median is another way to avoid extreme genes from swaying the value too much in One Direction genes with huge differences in expression have no more influence on the median than genes with minor differences since genes with huge differences will most likely be rare the effect is to give more influence to moderate differences and housekeeping genes okay now we're ready for step six convert the medians to normal numbers to get the final scaling factors for each sample again these are log values so they are exponents in this case exponents for E to calculate the scaling factor for each sample we raise e to the median value for each sample awesome we have scaling factors for the three samples now all we do is divide the original read counts by them that leads us to step seven divide the original read counts by the scaling factors here's our table that lists the original read counts note Gene one is part of the original read counts we didn't use Gene one when we calculated the scaling factors because it has a zero for read counts in Sample number one however we still need to scale the read counts for the other samples and here we have a table of scaled read counts I've rounded to the nearest read just to make this table easy to look at we can see that the read counts for sample number one were scaled up and that the read counts for sample number three were scaled down here's a summary of De 2's Library size scaling Factor logs eliminate all genes that are only transcribed in one sample type liver versus spleen they also help smooth over outlier read counts via the geometric mean the median further downplays genes that se soak up a lot of reads putting more emphasis on moderately expressed genes the idea behind using logs and the median is to hopefully Focus the scaling Factor on just the housekeeping genes the genes that are transcribed at the same levels in all of the samples you're looking at hooray we've made it to the end we now know how de 2 normalizes the read counts in each Library tune in next time when we talk about how Edge R normalizes the read counts until then Quest on