Calculate transcript abundances

Stringtie

Here we’ll import the Stringtie data:

# For each strain
for (strain in c("3d7","hb3","it")) {
  # create a tibble with the right columns
  assign(x     = paste0("x",strain,"_tpm"),
         value = tibble::tibble(transcript_id=character(),
                      cov=double(),
                      FPKM=double(),
                      TPM=double(),
                      tp=integer(),
                      strain=character()))
  # for each timepoint
  for (tp in seq(1,7)) {
    # generate file path
    f <- paste0("../output/transcript_abundance/stringtie/",strain,"_tp",tp,"_rf.gtf")
    # import the file and convert to the right format
    df <- tibble::as_tibble(import.gff(f)) %>% 
            dplyr::filter(type=="transcript") %>% 
            dplyr::select(transcript_id,cov,FPKM,TPM) %>% 
            dplyr::mutate(cov=as.double(cov),FPKM=as.double(FPKM),TPM=as.double(TPM),tp=tp,strain=strain) %>%
            dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"rna_","")) %>%
            dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"-1$",""))
    # bind this data to the tibble we created earlier for each strain
    assign(x     = paste0("x",strain,"_tpm"),
           value = dplyr::bind_rows(eval(parse(text=paste0("x",strain,"_tpm"))),df))
  }
}
# clean up
rm(df)

And let’s save these to output:

# Sense and antisense counts in long format
saveRDS(object=x3d7_tpm,file="../output/transcript_abundance/stringtie_3d7_abund.rds")
saveRDS(object=xhb3_tpm,file="../output/transcript_abundance/stringtie_hb3_abund.rds")
saveRDS(object=xit_tpm,file="../output/transcript_abundance/stringtie_it_abund.rds")

HTSeq

In order to properly calculate TPMs, we need to calculate the effective transcript lengths. To do this we need the mean or median fragment length distribution for each sample. We can generate this data using picard:

for file in $(find data/bam/mapped_to_3d7/ -name *.bam); do 
  picard CollectInsertSizeMetrics INPUT=$file HISTOGRAM_FILE=${file%.*}_insert_size.hist OUTPUT=${file%.*}_insert_size.out; 
done

Now we can read each file in and record the mean fragment length:

# create a tibble with the right columns
inserts <- tibble(sample=character(),mean_insert_size=double(),median_insert_size=double(),tp=integer())
# for each strain
for (strain in c("3d7","hb3","it")) {
  for (tp in seq(1,7)) {
    # generate the file path
    f <- paste0("../data/bam/mapped_to_3d7/",strain,".3d7_v3_chr.tp",tp,"_insert_size.out")
    # extract mean and median insert sizes
    mean_is <- read_tsv(f,skip=6,n_max = 1)$MEAN_INSERT_SIZE
    median_is <- read_tsv(f,skip=6,n_max = 1)$MEDIAN_INSERT_SIZE
    # bind this row to the tibble we created earlier
    inserts <- dplyr::bind_rows(inserts,tibble(sample=paste0("x",strain,"_tp",tp),mean_insert_size=mean_is,median_insert_size=median_is,tp=tp))
  }
}

It turns out, however, that using effective transcript lengths in Plasmodium is problematic…Some transcripts end up with a negative length if we use the formula provided for us here. So we’ll simply use full transcript lengths instead of “correcting” for that bias.

Import the HTSeq data:

# import our calculate transcript lengths
transcript_lengths <- read_tsv("../output/transcript_abundance/transcript_lengths.tsv")

# for each strain
for (strain in c("3d7","hb3","it")) {
  # create a tibble with the right colums
  assign(x     = paste0("x",strain,"_counts"),
         value = tibble(transcript_id=character(),
                        sense_counts=double(),
                        antisense_counts=double()))
  # for each time point
  for (tp in seq(1,7)) {
    # generate the file names
    s <- paste0("../output/transcript_abundance/htseq/",strain,"_tp",tp,"_reverse.txt")
    a <- paste0("../output/transcript_abundance/htseq/",strain,"_tp",tp,"_yes.txt")
    # and import sense and antisense counts and process their raw values
    st <- read_tsv(s,comment="__",col_names=c("transcript_id","sense_counts")) %>%
      dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"rna_","")) %>%
      dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"-1$",""))
    at <- read_tsv(a,comment="__",col_names=c("transcript_id","antisense_counts")) %>%       
      dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"rna_","")) %>%
      dplyr::mutate(transcript_id=stringr::str_replace(transcript_id,"-1$",""))
    # join the sense and antisense counts and give it an extra column
    df <- inner_join(st,at) %>% mutate(tp=tp)
    # bind the rows to the strain specific tibble generated above
    assign(x     = paste0("x",strain,"_counts"),
           value = dplyr::bind_rows(eval(parse(text=paste0("x",strain,"_counts"))),df))
  }
  # join this tibble with our transcript lengths for convenient TPM calculations
   assign(x     = paste0("x",strain,"_counts"),
           value = dplyr::inner_join(eval(parse(text=paste0("x",strain,"_counts"))),transcript_lengths))
}
# clean up
rm(s,a,st,at,df)

Calculate TPM values:

# Function for calculating TPMs from the count tibbles
# generated above
# 
calc_tpms <- function(counts) {
  # caluclate the normalizing factor for each sample
  # normFactor = sum of all length normalized transcript counts
  normFactors <- counts %>% 
  mutate(normCount=sense_counts/length) %>% 
  group_by(tp) %>% 
  summarise(normFactor=1/sum(normCount))
  
  # calculate the TPMs by joining it together with normFactors
  # and multiplying the length corrected counts by these factors
  # and by one million
  abund <- inner_join(normFactors,counts) %>%
  mutate(sense_tpm=(sense_counts/length)*normFactor*1e6,antisense_tpm=(antisense_counts/length)*normFactor*1e6) %>%
  dplyr::select(transcript_id,sense_counts,antisense_counts,sense_tpm,antisense_tpm,tp)
  
  return(abund)
}

# calculate the TPMs for each strain
x3d7_abund <- calc_tpms(x3d7_counts)
xhb3_abund <- calc_tpms(xhb3_counts)
xit_abund  <- calc_tpms(xit_counts)

Let’s just make sure they all add up to a million as they should:

# sum up TPM values for each time point
x3d7_abund %>% group_by(tp) %>% summarise(total=sum(sense_tpm))

# A tibble: 7 × 2
     tp total
  <int> <dbl>
1     1 1e+06
2     2 1e+06
3     3 1e+06
4     4 1e+06
5     5 1e+06
6     6 1e+06
7     7 1e+06

xhb3_abund %>% group_by(tp) %>% summarise(total=sum(sense_tpm))

# A tibble: 7 × 2
     tp total
  <int> <dbl>
1     1 1e+06
2     2 1e+06
3     3 1e+06
4     4 1e+06
5     5 1e+06
6     6 1e+06
7     7 1e+06

xit_abund %>% group_by(tp) %>% summarise(total=sum(sense_tpm))

# A tibble: 7 × 2
     tp total
  <int> <dbl>
1     1 1e+06
2     2 1e+06
3     3 1e+06
4     4 1e+06
5     5 1e+06
6     6 1e+06
7     7 1e+06

Write these abundance values to output:

# Sense TPM tables
select(x3d7_abund,transcript_id,sense_tpm,tp) %>% spread(tp,sense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_3d7_sense_tpms.tsv")
select(xhb3_abund,transcript_id,sense_tpm,tp) %>% spread(tp,sense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_hb3_sense_tpms.tsv")
select(xit_abund,transcript_id,sense_tpm,tp) %>% spread(tp,sense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_it_sense_tpms.tsv")

# Antisense TPM tables
select(x3d7_abund,transcript_id,antisense_tpm,tp) %>% spread(tp,antisense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_3d7_antisense_tpms.tsv")
select(xhb3_abund,transcript_id,antisense_tpm,tp) %>% spread(tp,antisense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_hb3_antisense_tpms.tsv")
select(xit_abund,transcript_id,antisense_tpm,tp) %>% spread(tp,antisense_tpm) %>%
  write_tsv(path="../output/transcript_abundance/htseq_it_antisense_tpms.tsv")

# Sense and antisense counts in long format
saveRDS(object=x3d7_abund,file="../output/transcript_abundance/htseq_3d7_abund.rds")
saveRDS(object=xhb3_abund,file="../output/transcript_abundance/htseq_hb3_abund.rds")
saveRDS(object=xit_abund,file="../output/transcript_abundance/htseq_it_abund.rds")

That’s it! Now we can use these files for future analyses.