Spen_RNAseq_comparison.Rmd

---
title: "S. pennellii RNAseq simulations"
author: "Julin Maloof"
date: "July 26, 2016"
output:
  revealjs::revealjs_presentation:
    reveal_options:
      width: 1200
      height: 800
---

```{r setup, include=FALSE}
verbose=FALSE
library(knitr)
opts_chunk$set(cache.lazy=FALSE)
opts_chunk$set(cache=TRUE)
opts_chunk$set(warnings=FALSE)
opts_chunk$set(messages=FALSE)
```

```{r load_libraries, cache=FALSE, include=verbose}
library(polyester)
library(Biostrings)
library(ggplot2)
library(Rsamtools)
library(rtracklayer)
library(Rsubread)
library(R.utils)
library(plyr)
library(magrittr)
library(reshape2)
setwd("~/Documents/Lab Notebook support/2015/rnaseq_simulations")
source("Rnaseq_sim_helper_functions.R")
```

# Introduction

## Overall Question
<!-- .slide: style="text-align: left;"> --> 
When mapping RNAseq reads from a heterologous species to a reference genome, what is the best way to perform the mapping?

For example:

> - We have a tomato reference genome
> - We want to do RNAseq on a variety of wild relatives of tomato for which we have no genome

## Strategy: simulation

> * Use _S. pennellii_ and _S. lycopersium_.  Both are sequenced.  
> * Simulate RNASeq reads from _S. pennellii_ and map to _S. lycopersicum_.
> * Compare actual (simulated) counts to the number of counts that map.

> * Can compare different aligners:
> * Splicing aligners that map RNAseq to genome
    * [stampy](http://www.well.ox.ac.uk/project-stampy)
    * [subread](http://bioinf.wehi.edu.au/subread/)
    * [STAR](https://github.com/alexdobin/STAR)
    * [bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/index.shtml)
> * Non-splicing aligners that map to a cDNA reference
    * [BWA](http://bio-bwa.sourceforge.net/)
    * [Kallisto](https://pachterlab.github.io/kallisto/)

> * And different parameters with those aligners.

# Simulate Reads

## Use Polyester

Use the [Polyester package](http://bioconductor.org/packages/release/bioc/vignettes/polyester/inst/doc/polyester.html)  to simulate RNAseq reads.  Also see the [Polyester paper.](http://bioinformatics.oxfordjournals.org/content/31/17/2778)

## Simluate reads (simple version)
Load Spen and Slyc transcripts  
Spen CDS fasta can be obtained from [SGN](ftp://ftp.solgenomics.net/genomes/Solanum_pennellii/annotations/)
slyc fasta is also from [SGN](ftp://ftp.sgn.cornell.edu/genomes/Solanum_lycopersicum/annotation/ITAG2.4_release/)

```{r load_transcripts, eval=TRUE, include=verbose}
pen.transcripts <-  readDNAStringSet("~/Sequences/ref_genomes/S.pen/Spenn-v2-cds-annot.fa")
head(pen.transcripts)
names(pen.transcripts[1:10])
lyc.transcripts <- readDNAStringSet("~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/ITAG2.4_cds.fasta")
head(lyc.transcripts)
names(lyc.transcripts)[1:10]
```

## Reformat names to get rid of description and add Slyc reference
```{r reformat_gene_names, eval=TRUE}
orthologs <- read.delim("Ortho_v2.txt",header=FALSE) #From Tony Bolger
head(orthologs)
names(orthologs) <- c("PEN","LYC")
orthologs$PEN <- sub("t","",orthologs$PEN)
orthologs$PEN.LYC <- paste(orthologs$PEN,orthologs$LYC,sep="_")
names(pen.transcripts) <- substr(names(pen.transcripts),1,16)
names(pen.transcripts) <- ifelse(!is.na(match(names(pen.transcripts),orthologs$PEN)),
                                 orthologs$PEN.LYC[match(names(pen.transcripts),orthologs$PEN)],
                                 names(pen.transcripts))
pen.transcripts

names(lyc.transcripts) <- substr(names(lyc.transcripts),1,18)
names(lyc.transcripts) <- ifelse(!is.na(match(names(lyc.transcripts),orthologs$LYC)),
                                 orthologs$PEN.LYC[match(names(lyc.transcripts),orthologs$LYC)],
                                 names(lyc.transcripts))
pen.gene.length <- as.integer(nchar(pen.transcripts))
names(pen.gene.length) <- names(pen.transcripts)
lyc.gene.length <- as.integer(nchar(lyc.transcripts))
names(lyc.gene.length) <- names(lyc.transcripts)
```

## Limit ourselves to genes where orthology is clear:
```{r reduce_data, eval=TRUE, include=verbose}
pen.transcripts.small <- pen.transcripts[nchar(names(pen.transcripts))>16]
lyc.transcripts.small <- lyc.transcripts[nchar(names(lyc.transcripts)) > 18]

if(! file.exists("~/Sequences/ref_genomes/S.pen/Spenn-v2-cds-annot_orthlogs_only.fa")) writeXStringSet(pen.transcripts.small,file="~/Sequences/ref_genomes/S.pen/Spenn-v2-cds-annot_orthlogs_only.fa")

if(! file.exists("~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/ITAG2.4_cds.pen.orthologs.only.fasta")) writeXStringSet(lyc.transcripts.small,file="~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/ITAG2.4_cds.pen.orthologs.only.fasta")
```

## Simulate PEN RNAseq reads

Download [fasta to fastq converter](https://code.google.com/archive/p/fasta-to-fastq/downloads)

```{r simulate_pen_3, eval=FALSE}
#4.25 hrs for 6 * 10M reads
system.time(
  simulate_experiment(fasta="~/Sequences/ref_genomes/S.pen/Spenn-v2-cds-annot_orthlogs_only.fa",
                      outdir="spen_sim_3",
                      num_reps=6,
                      readlen=50,
                      #line below gives ~ 10M reads, scaled by transcript length and following a rexp distribution
                      reads_per_transcript = round(rexp(n=length(pen.transcripts.small)) *
                                                     width(pen.transcripts.small) * 
                                                     1e7/ sum(width(pen.transcripts.small) ))+1,
                      paired=FALSE,
                      error_model="illumina5",
                      bias="rnaf"
  )
)
system("for f in `ls spen_sim_3/sample_*`
do
echo $f
fnew=`echo $f | sed s/.fasta/_renamed.fasta/`
fnewq=`echo $f | sed s/.fasta/_renamed.fastq/`
sed 's!/!;!' $f > $fnew #allows original gene_names to be kept
rm $f
cat  $fnew | perl ~/bin/fasta_to_fastq.pl  - | gzip > $fnewq.gz
gzip $fnew
done")

```

# Map

## setup for mapping

Default parameters for everything but kallisto, which uses kmer=19
```{r, eval=FALSE, include=verbose}
library(snowfall)
files <- dir(path="./spen_sim_3",pattern="fasta.gz")
filesq <- dir(path="./spen_sim_3",pattern="fastq.gz")
sfInit(parallel = T, cpus = min(length(files),12))
sfSource("Rnaseq_sim_helper_functions.R")
sfLibrary(Rsubread)
sfLibrary(R.utils)
```

## Star
```{r, eval=FALSE}
runStar.time <- sfSapply(files,function(f) {
  system.time(RunStar(fasta=f,
                      prefix = paste("STAR",substr(f,1,9),sep="_"),
                      param.file = "../STAR.params.whitney.1", 
                      dir="spen_sim_3",
                      n=2))
})
```

## Kallisto

```{r, eval=FALSE}
runKallisto.time <- sfSapply(filesq,function(f) {
  system.time(RunKallisto(fastq=f,
                          index = "ITAG2.4._cds.pen.orthlogsonly.19.kai",
                          dir = "spen_sim_3",
                          prefix=paste("kallisto",substr(f,1,9),sep="_"),
                          threads=2))
})
```

## Subread

```{r, eval=FALSE}
runSubread.time <- sfSapply(filesq,function(f) {
  system.time(align("Slyc_Chromo2.5_subread_index",
                    file.path("spen_sim_3",f), 
                    output_file=paste("spen_sim_3/subread_default_",substr(f,1,9),".bam",sep=""),
                    nthreads = 2))
})
```

## Bowtie2

```{r, eval=FALSE}
runBowtie2.time <- sfSapply(filesq, function(f) {
  system.time(RunBowtie2(fastq = f,
                         index = "~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/S_lycopersicum_chromosomes.2.50.fa",
                         dir = "spen_sim_3/",
                         prefix=paste("bowtie2",substr(f,1,9),sep="_"),
                         threads=2))
})
```

## BWA

```{r, eval=FALSE}
runBWA.time <- sfSapply(files, function(f) {
  system.time(RunBWAaln(f,
                        index = "~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/ITAG2.4_cds.pen.orthologs.only.fasta",
                        dir = "spen_sim_3", 
                        prefix=paste("BWA",substr(f,1,9),sep="_"),
                        threads=2))
})
```

## Stampy

```{r, eval=FALSE}
runStampy.time <- sfSapply(files,function(f) {
  system.time(RunStampy(f,
                        index = "SL2.5",
                        dir = "spen_sim_3", 
                        prefix=paste("stampy",substr(f,1,9),sep="_"),
                        threads=2))
})
```

## save it

```{r, eval=FALSE}
save.image("runmapstats.Rdata")
```

#get mapped reads

## Star counts

```{r, eval=FALSE}
star.bams <- dir("spen_sim_3",pattern="STAR_.+bam$")
star.counts.time <- system.time(star.mapped.counts <- sfSapply(star.bams, function(bam) {
  GetMappedCounts(bam = bam, dir = "spen_sim_3", type="featureCounts")
}))
```

## Kallisto counts

```{r, eval=FALSE}
kallisto.dirs <- dir("spen_sim_3",pattern="kallisto_sample")
kallisto.counts.time <- system.time(kallisto.mapped.counts <- sfSapply(kallisto.dirs, function(d) {
  GetMappedCounts(dir = file.path("spen_sim_3",d), type="kallisto-transcripts")
}))
```

## Subread counts

```{r, eval=FALSE}
subread.bams <- dir("spen_sim_3",pattern="subread_.+bam$")
subread.counts.time <- system.time(subread.mapped.counts <- sfSapply(subread.bams, function(bam) {
  GetMappedCounts(bam = bam, dir = "spen_sim_3", type="featureCounts")
}))
```

## Bowtie2 counts

```{r, eval=FALSE}
bowtie2.bams <- dir("spen_sim_3",pattern="bowtie2_sample_.+bam$")
bowtie2.counts.time <- system.time(bowtie.mapped.counts <- sfSapply(bowtie2.bams, function(bam) {
  GetMappedCounts(bam = bam, dir = "spen_sim_3", type="featureCounts")
}))
```

## BWA counts

```{r, eval=FALSE}
bwa.bams <- dir("spen_sim_3",pattern="BWA_sample_.+bam$")
bwa.counts.time <- system.time(bwa.mapped.counts <- sfSapply(bwa.bams, function(bam) {
  GetMappedCounts(dir="spen_sim_3/",type = "transcripts",bam = bam)
}))
```

## stampy counts

```{r, eval=FALSE}
stampy.bams <- dir("spen_sim_3",pattern="stampy_sample_")
stampy.counts.time <- system.time(stampy.mapped.counts <- sfSapply(stampy.bams, function(bam) {
  GetMappedCounts(bam = bam, dir = "spen_sim_3", type="featureCounts")
}))

save.image("counts_and_times.Rdata")
```

# plot comparisons

## compare run time

```{r, echo=FALSE}
program.order <- c("Bowtie2","Stampy","Star","Subread","BWA","Kallisto")
load("counts_and_times.Rdata")
run.times <- ls(pattern="run.+time$")
time.data <- sapply(run.times,function(run) {
  cbind(mean(get(run)["elapsed",]),
        sem(get(run)["elapsed",]))}) %>%
  t() %>% 
  as.data.frame()
colnames(time.data) <- c("mean","sem")

time.data$program <- sub("run","",rownames(time.data)) %>% 
  sub(".time","",.) %>%
  factor(., levels=program.order)

time.pl <- ggplot(time.data,aes(x=program,y=mean,ymin=mean-sem,ymax=mean+sem))
time.pl <- time.pl + geom_bar(stat="identity",fill="skyblue")
time.pl <- time.pl + geom_errorbar(width=.5)
time.pl + ggtitle("Average run time, 10M reads, 2 CPUs") + ylab("time (sec)")  +  theme(axis.text.x  = element_text(size=16),axis.text.y=element_text(size=16))
```

## compare % reads mapped

```{r get_known_counts,include=verbose, eval=FALSE}
pen.known.counts.list <- sfLapply(files, function(f) {
  tmp <- read.table(
    pipe(paste("./CountGene.py",file.path("spen_sim_3",f))),
    col.names=c("id","count"),stringsAsFactors = FALSE)
  tmp$lyc.id <- substr(tmp$id,18,33)
  tmp })

save(pen.known.counts.list,file="pen.known.counts.Rdata")
```

```{r known.counts.transform, include=verbose}
load("pen.known.counts.Rdata")

names(pen.known.counts.list) <- substr(files,1,9)

pen.known.counts.m <- ldply(pen.known.counts.list,.id="sample")

pen.known.counts <- dcast(pen.known.counts.m, id + lyc.id ~ sample, value.var = "count")

```

```{r get_reads, include=verbose}
bwa.counts.matrix <- apply(bwa.mapped.counts, 2, function(x) x$count)
row.names(bwa.counts.matrix) <- bwa.mapped.counts[,1]$ID

bowtie.counts.matrix <- apply(bowtie.mapped.counts,2, function(x) x$count)
row.names(bowtie.counts.matrix) <- bowtie.mapped.counts[,1]$ID
bowtie.counts.matrix <- bowtie.counts.matrix[grepl("Solyc",row.names(bowtie.counts.matrix)),] # get rid of rows with "stats" data

kallisto.counts.matrix <- apply(kallisto.mapped.counts,2, function(x) x$count)
row.names(kallisto.counts.matrix) <- kallisto.mapped.counts[,1]$ID
kallisto.counts.matrix <- kallisto.counts.matrix[grepl("Solyc",row.names(kallisto.counts.matrix)),] # get rid of rows with "stats" data

stampy.counts.matrix <- apply(stampy.mapped.counts,2, function(x) x$count)
row.names(stampy.counts.matrix) <- stampy.mapped.counts[,1]$ID
stampy.counts.matrix <- stampy.counts.matrix[grepl("Solyc",row.names(stampy.counts.matrix)),] # get rid of rows with "stats" data

star.counts.matrix <- apply(star.mapped.counts,2, function(x) x$count)
row.names(star.counts.matrix) <- star.mapped.counts[,1]$ID
star.counts.matrix <- star.counts.matrix[grepl("Solyc",row.names(star.counts.matrix)),] # get rid of rows with "stats" data

subread.counts.matrix <- apply(subread.mapped.counts,2, function(x) x$count)
row.names(subread.counts.matrix) <- subread.mapped.counts[,1]$ID
subread.counts.matrix <- subread.counts.matrix[grepl("Solyc",row.names(subread.counts.matrix)),] # get rid of rows with "stats" data
```

```{r compare_percent_mapped, echo=FALSE}
total.reads <- tapply(pen.known.counts.m$count,pen.known.counts.m$sample,sum)

counts.objects <- ls(pattern="counts\\.matrix$")

percent.mapped.data <- sapply(counts.objects, function(x) {
  cbind(mean(colSums(get(x)) / total.reads * 100),
        sem(colSums(get(x)) / total.reads * 100))}) %>%
  t() %>%
  as.data.frame()

colnames(percent.mapped.data) <- c("mean","sem")

percent.mapped.data$program <- factor(sub("\\.counts\\.matrix","",row.names(percent.mapped.data)),
                                      levels=c("bowtie","stampy","star","subread","bwa","kallisto"),
                                          labels=c("Bowtie2","Stampy","Star","Subread","BWA","Kallisto"))

percent.map.pl <- ggplot(percent.mapped.data,aes(x=program,y=mean,ymin=mean-sem,ymax=mean+sem))
percent.map.pl <- percent.map.pl + geom_bar(stat="identity",fill="skyblue")
percent.map.pl <- percent.map.pl + geom_errorbar(width=.5)
percent.map.pl + ggtitle("Perecent Mapped Reads") + ylab("Percent") + coord_cartesian(ylim=c(90,100)) +  theme(axis.text.x  = element_text(size=16),axis.text.y=element_text(size=16))

```


## compare correlations all genes (Pearson)

```{r, echo=FALSE,results='hide'}
samples <- paste("sample_0",as.character(1:6),sep="")

cor.data <- t(sapply(ls(pattern="counts.matrix$"),function(m) {
  tmp.m <- get(m)

tmp.df <- data.frame(
  ID=row.names(tmp.m),
  tmp.m)

cors <- as.data.frame(t(sapply(samples, function(s) {
  CompareCounts(known.counts=subset(pen.known.counts,select=c("lyc.id",grep(s,colnames(pen.known.counts),value=TRUE))),
                mapped.counts=subset(tmp.df,select=c("ID",grep(s,colnames(tmp.df),value=TRUE))),
                plot = FALSE,return.merged.table = FALSE)
})))
c(Pearson.mean=mean(cors$Pearson),
  Pearson.sem=sem(cors$Pearson),
  Spearman.mean=mean(cors$Spearman),
  Spearman.sem=sem(cors$Spearman))
}
))
cor.data <- as.data.frame(cor.data)

cor.data$program <- factor(sub("\\.counts\\.matrix","",row.names(cor.data)),
                                      levels=c("bowtie","stampy","star","subread","bwa","kallisto"),
                                          labels=c("Bowtie2","Stampy","Star","Subread","BWA","Kallisto"))

cor.pl <- ggplot(data=cor.data,aes(x=program,y=Pearson.mean,ymin=Pearson.mean+Pearson.sem,ymax=Pearson.mean+Pearson.sem))
cor.pl <- cor.pl + geom_bar(stat="identity",fill="skyblue")
cor.pl <- cor.pl + geom_errorbar(width=0.5)
cor.pl + ylab("mean Pearson correlation") + coord_cartesian(ylim=c(0.9,1)) + ggtitle("Pearson Correlation")  +  theme(axis.text.x  = element_text(size=16),axis.text.y=element_text(size=16))
```

## compare correlations all genes (Spearman)

```{r, echo=FALSE}
cor.pl <- ggplot(data=cor.data,aes(x=program,y=Spearman.mean,ymin=Spearman.mean+Spearman.sem,ymax=Spearman.mean+Spearman.sem))
cor.pl <- cor.pl + geom_bar(stat="identity",fill="skyblue")
cor.pl <- cor.pl + geom_errorbar(width=0.5)
cor.pl + ylab("mean Spearman correlation") + coord_cartesian(ylim=c(0.9,1)) + ggtitle("Spearman Correlation")  +  theme(axis.text.x  = element_text(size=16),axis.text.y=element_text(size=16))
```



## correlation plots all genes
```{r, echo=FALSE}
op <- par(mfrow=c(3,2))

sapply(ls(pattern="counts.matrix$")[c(1,4,5,6,2,3)],function(m) {
  tmp.m <- get(m)
  tmp.df <- data.frame(
    ID=row.names(tmp.m),
    tmp.m)

  CompareCounts(known.counts=subset(pen.known.counts,select=c("lyc.id",grep("sample_01",colnames(pen.known.counts),value=TRUE))),
                mapped.counts=subset(tmp.df,select=c("ID",grep("sample_01",colnames(tmp.df),value=TRUE))),
                plot = TRUE, chrom.separate=FALSE,return.merged.table = FALSE, correlation = FALSE, title = sub("\\.counts\\.matrix","",m))
})
par(op)
```


## compare correlations "bad" genes removed

```{r, echo=FALSE}
#first find genes with overlapping annotation
lyc.transcript.gr <- import("~/Sequences/ref_genomes/tomato/ITAG2.4_Chromo2.5/ITAG2.4_gene_models.gff3",feature.type="mRNA")

lyc.lyc.overlaps <- findOverlaps(lyc.transcript.gr,lyc.transcript.gr) #all overlaps including self on self
lyc.lyc.overlaps <- lyc.lyc.overlaps[queryHits(lyc.lyc.overlaps)!=subjectHits(lyc.lyc.overlaps),] #limit to non-self overlaps
lyc.lyc.overlaps <- substr(unique(mcols(lyc.transcript.gr)[queryHits(lyc.lyc.overlaps),"Name"]),1,16)

pen.transcript.gr <- PSL2Granges(file="penCDS_to_lycGENME.psl")
pen.pen.overlaps <- findOverlaps(pen.transcript.gr,pen.transcript.gr) #all overlaps including self on self

#need additional step to limit to non-self overlaps because BLAT returns multiple hits in the same region.  Make sure that the overlap is with the same gene
pen.pen.overlaps <- pen.pen.overlaps[names(pen.transcript.gr)[queryHits(pen.pen.overlaps)]!=names(pen.transcript.gr)[subjectHits(pen.pen.overlaps)]] 

pen.pen.overlaps <- unique(mcols(pen.transcript.gr)[queryHits(pen.pen.overlaps),"lyc.ID"]) 
```

```{r, echo=FALSE,results='hide'}
pen.known.counts.no.overlaps <- pen.known.counts[! pen.known.counts$lyc.id %in% c(pen.pen.overlaps,lyc.lyc.overlaps),]

samples <- paste("sample_0",as.character(1:6),sep="")

cor.data <- t(sapply(ls(pattern="counts.matrix$"),function(m) {
  tmp.m <- get(m)

tmp.df <- data.frame(
  ID=row.names(tmp.m),
  tmp.m)
tmp.df <- tmp.df[! tmp.df$ID %in% c(pen.pen.overlaps,lyc.lyc.overlaps),]

cors <- as.data.frame(t(sapply(samples, function(s) {
  CompareCounts(known.counts=subset(pen.known.counts.no.overlaps,select=c("lyc.id",grep(s,colnames(pen.known.counts.no.overlaps),value=TRUE))),
                mapped.counts=subset(tmp.df,select=c("ID",grep(s,colnames(tmp.df),value=TRUE))),
                plot = FALSE,return.merged.table = FALSE)
})))
c(Pearson.mean=mean(cors$Pearson),
  Pearson.sem=sem(cors$Pearson),
  Spearman.mean=mean(cors$Spearman),
  Spearman.sem=sem(cors$Spearman))
}
))
cor.data <- as.data.frame(cor.data)

cor.data$program <- factor(sub("\\.counts\\.matrix","",row.names(cor.data)),
                                      levels=c("bowtie","stampy","star","subread","bwa","kallisto"),
                                          labels=c("Bowtie2","Stampy","Star","Subread","BWA","Kallisto"))

cor.pl <- ggplot(data=cor.data,aes(x=program,y=Spearman.mean,ymin=Spearman.mean+Spearman.sem,ymax=Spearman.mean+Spearman.sem))
cor.pl <- cor.pl + geom_bar(stat="identity",fill="skyblue")
cor.pl <- cor.pl + geom_errorbar(width=0.5)
cor.pl + ylab("mean Spearman correlation") + coord_cartesian(ylim=c(0.9,1)) + ggtitle("Spearman Correlation")  +  theme(axis.text.x  = element_text(size=16),axis.text.y=element_text(size=16))
```

## correlation plots "bad" genes removed
```{r, echo=FALSE}
op <- par(mfrow=c(3,2))

sapply(ls(pattern="counts.matrix$")[c(1,4,5,6,2,3)],function(m) {
  tmp.m <- get(m)
  tmp.df <- data.frame(
    ID=row.names(tmp.m),
    tmp.m)
  tmp.df <- tmp.df[! tmp.df$ID %in% c(pen.pen.overlaps,lyc.lyc.overlaps),]
  
  CompareCounts(known.counts=subset(pen.known.counts.no.overlaps,select=c("lyc.id",grep("sample_01",colnames(pen.known.counts.no.overlaps),value=TRUE))),
                mapped.counts=subset(tmp.df,select=c("ID",grep("sample_01",colnames(tmp.df),value=TRUE))),
                plot = TRUE, chrom.separate=FALSE,return.merged.table = FALSE, correlation = FALSE, title = sub("\\.counts\\.matrix","",m))
})
par(op)
```


# ToDos

> * Check for correct splicing across exons (IMPORTANT for SNPs)
> * How many reads are mapped to the wrong gene?
> * How many genes have a counts ratio < 0.9 or > 1.1
> * Add more SNPs/inDels
> * Check for correct splicing across exons (IMPORTANT for SNPs)