Analyzing compositional count data from Maize TIR experiment

Load data from counting scripts

library(compositions)
# read library totals
tot<- read.csv("~/Documents/collaboratorProjects/chip_hunter_maize_TIR/data/totals.txt", sep ="\t", header=F)

Option 1: Close the composition using the total reads

These steps perform the closure operation telling it the total number of reads, then cast the data as a count composition dataframe that sums to 1. Finally we calculate a centered log ratio transformation on the data.

compclosed <- clo(X=a2t,total=tot[,1], detectionlimit= 1)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

# create a ccount composition dataset
a2c<-ccomp(compclosed, total=1)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

#calulate a centered log ratio transform 
a2d <- clr(a2c)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

Heatmap

A heat map of all data shows some general trends (forward and reverse probes tend to group together)

# plot a heat map of the results
library(gplots)

Attaching package: ‘gplots’

The following object is masked from ‘package:stats’:

    lowess

library(viridis)

Loading required package: viridisLite

heatmap.2(a2d, col=viridis, scale="none",trace="none")

Analysis of known homozygotes and heterozygotes

Next we attempt to plot the forward and reverse data for two probes where we know Homozygous samples are in 1-15 and heterozygous are in 18-30. This data is the clr values for the dataset closed on all reads.

plot(a2d[,85])

plot(a2d[,86])

fr = a2d[,85] + a2d[,86]
plot(fr)# Plot mean and SD of Homo

mean(fr[1:15])

[1] 0.7546468

msd(fr[1:15])

[1] 1.601614

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] 0.401688

msd(fr[18:30])

[1] 1.581361

The data is not very distinguishable. Other normalization approaches are needed.

Option 2: Close the composition using the total mapped reads

compclosed2 <- clo(X=a2t, detectionlimit= 1)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

# create a ccount composition dataset
a2c2<-ccomp(compclosed2, total=1)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

#calulate a centered log ratio transform 
a2d2 <- clr(a2c2)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

Analysis of known homozygotes and heterozygotes

Next we attempt to plot the forward and reverse data for two probes where we know Homozygous samples are in 1-15 and heterozygous are in 18-30.

This data is the count composition value for each probe

plot(a2c2[,85])

plot(a2c2[,86])

fr = a2c2[,85] + a2c2[,86]
plot(fr)

# Plot mean and SD of Homo
mean(fr[1:15])

[1] 0.03628429

msd(fr[1:15])

[1] 0.02337816

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] 0.03999595

msd(fr[18:30])

[1] 0.02945496

This data is the clr transformation of the count composition value for each probe

# plot the clr data for probes "Probe087-mu1017716-F", "Probe088-mu1017716-R" where we know half are heterozygous and half are homozugous
plot(a2d2[,85])

plot(a2d2[,86])

fr = a2d2[,85] + a2d2[,86]
plot(fr)

# Plot mean and SD of Homo
mean(fr[1:15])

[1] 0.7546468

msd(fr[1:15])

[1] 1.601614

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] 0.401688

msd(fr[18:30])

[1] 1.581361

Neither approach is easily separable.

Option 3: treat the unmapped reads as a component of the composition

cs <-colSums(a2)
nonmapped <- tot$reads-cs
a3 <- rbind(a2,nonmapped)
rownames(a3)[145]<-"nonmapped"
a3t=t(a3)
a3c<-ccomp(a3t, total=1)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

a3d <- clr(a3c)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

Heatmap

A heat map of all data shows some general trends (forward and reverse probes tend to group together)

# plot a heat map of the results
library(gplots)
library(viridis)
heatmap.2(a3d, col=viridis, scale="none",trace="none")

Analysis of known homozygotes and heterozygotes

plot(a3c[,85])

plot(a3c[,86])

fr = a3c[,85] + a3c[,86]
plot(fr)

# Plot mean and SD of Homo
mean(fr[1:15])

[1] 4.432075e-05

msd(fr[1:15])

[1] 7.462433e-05

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] 3.889691e-05

msd(fr[18:30])

[1] 4.013475e-05

These results are closer being seperable

plot the clr of the compositional data for probes “Probe087-mu1017716-F”, “Probe088-mu1017716-R” where we know half are heterozygous and half are homozugous

plot(a3d[,85])

plot(a3d[,86])

fr = a3d[,85] + a3d[,86]
plot(fr)

# Plot mean and SD of Homo
mean(fr[1:15])

[1] 1.473779

msd(fr[1:15])

[1] 2.061553

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] 1.415554

msd(fr[18:30])

[1] 2.028346

the trend is less clear

Option 4: Additive log ratio

Use the additive log ratio with non-mapped reads as the common denominator

a3a <- alr(a3c,ivar=145)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

plot the compositional data for probes “Probe087-mu1017716-F”, “Probe088-mu1017716-R” where we know half are heterozygous and half are homozugous

plot(a3a[,85])

plot(a3a[,86])

fr = a3a[,85] + a3a[,86]
plot(fr)

# Plot mean and SD of Homo
mean(fr[1:15])

[1] -Inf

msd(fr[1:15])

[1] NaN

#Plot mean and SD of Hetero
mean(fr[18:30])

[1] -Inf

msd(fr[18:30])

[1] NaN

Again the means appear to be different but there is not a clear separation

Option 5: Balances

Balances are the geometric means of log ratios between sub-components of the data. We take the dataset where untapped is a component and we compare the log ratios of the forward probe to un-mapped. Again the data is not separable

a3acomp <- acomp(a3t)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

bal85 <- balance(X=acomp(a3t),expr= ~ (Probe087-mu1017716-F) / nonmapped)

Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.Calling 'structure(NULL, *)' is deprecated, as NULL cannot have attributes.
  Consider 'structure(list(), *)' instead.

plot(bal85[,1])

Option 6: simple normalization

For comparison here is simple normalization

simplenorm <- sweep(a2c, 2, tot[,1], "/")

STATS does not recycle exactly across MARGIN

plot(simplenorm[,85])

plot(simplenorm[,86])

Results

The best approach seems to be creating a compositional count dataset with the non mapped reads as a component. However, the difference is not enough to reliably classify the data.

Option 8: normalize by tircounts

For comparison here is simple normalization

tirnorm <- sweep(a2c, 2, tirlong, "/")

STATS does not recycle exactly across MARGIN

plot(tirnorm[,85])

plot(tirnorm[,86])

Enrichment efficiency

plot(log(rowSums(a2t)/tot[,1]))

table(a2t[,30] < 5)

FALSE  TRUE 
   24     6 

Enrichment efficiency using tircounts as the numerator

plot((tirlong/tot[,1]))

mean(log(tirlong/tot[,1]))

[1] -3.258941

sd(log(tirlong/tot[,1]))

[1] 0.7283162