Classifying putative mitochondrial contigs

Shaun Jackman 2013-12-12

Load libraries

## Loading required package: grid
## Loading required package: lattice
## Loading required package: ggplot2
## 
## Attaching package: 'useful'
## 
## The following object is masked from 'package:stats':
## 
##     plot.ts

Read the data

comp <- read.table('pg29organelle-scaffolds.fa.comp.tab', header=TRUE)
idxstats <- subset(read.table('pg29organelle-scaffolds_D08CWACXX_6_merged.sort.idxstats.tab',
    header=TRUE), tid != '*')
fatab <- read.table('pg29organelle-scaffolds.fa.tab', header=TRUE)
data <- cbind(fatab, comp, idxstats)

Perform calculations

data$numACGT <- with(data, numA + numC + numG + numT)
data$GCcontent <- with(data, (numC + numG) / numACGT)
data$read.density <- with(data, numMapped / numACGT)

# mean.coverage seems to perform better than read.density
data$read.density <- data$mean.coverage

Visualize the data in one dimension

size.threshold <- 2000

densityplot(~numACGT, data, subset = numACGT >= size.threshold,
    scales = list(x = list(log = TRUE)))

densityplot(~GCcontent, data, subset = numACGT >= size.threshold)

densityplot(~read.density, data, subset = numACGT >= size.threshold,
    scales = list(x = list(log = TRUE)))

Visualize the data in two dimensions

splom(~ cbind(read.density, GCcontent, numACGT), data,
    subset = numACGT >= size.threshold)

#xyplot(mean.coverage ~ read.density, data, subset = numACGT >= size.threshold)
xyplot(read.density ~ GCcontent, data, subset = numACGT >= size.threshold)

xyplot(read.density ~ numACGT, data, subset = numACGT >= size.threshold)

xyplot(GCcontent ~ numACGT, data, subset = numACGT >= size.threshold)

hexbinplot(read.density ~ GCcontent, data, subset = numACGT >= size.threshold)

hexbinplot(read.density ~ numACGT, data, subset = numACGT >= size.threshold)

hexbinplot(GCcontent ~ numACGT, data, subset = numACGT >= size.threshold)

xyplot(read.density ~ numACGT, data, subset = numACGT >= size.threshold,
    main = 'Read density vs. scaffold size of scaffolds larger than 5 kbp',
    xlab = 'Scaffold size (bp)',
    ylab = 'Read density (reads/bp)')

xyplot(500 * read.density ~ numACGT, data, subset = numACGT >= size.threshold,
    main = 'Coverage vs. scaffold size\nfor scaffolds larger than 5 kbp',
    xlab = 'Scaffold size (bp)',
    ylab = 'Read coverage')

Visualize the data in three dimensions

xyplot(read.density ~ GCcontent, data,
    subset = numACGT >= size.threshold,
    cex = log(subset(data, numACGT >= size.threshold)$numACGT / size.threshold, base=10))

Clustering with k-means

data.subset <- subset(data, numACGT >= size.threshold, c('name', 'read.density', 'GCcontent', 'numACGT'))
data.subset.scaled <- with(data.subset, data.frame(
    read.density = scale(log(read.density)),
    GCcontent = scale(GCcontent),
    numACGT = scale(log(numACGT))))
clusters <- kmeans(
    subset(data.subset.scaled, select = c('read.density', 'GCcontent', 'numACGT')),
    centers = 2)

plot(clusters, data = data.subset.scaled)

xyplot(read.density ~ GCcontent, data.subset.scaled, group = clusters$cluster,
    #cex = data.subset.scaled$numACGT - min(data.subset.scaled$numACGT),
    auto.key=T)

data.clustered <- cbind(data.subset, clusters['cluster'])

xyplot(read.density ~ GCcontent, data.clustered, group = cluster,
    #cex = log(data.clustered$numACGT / size.threshold, base=10),
    auto.key=T)

xyplot(read.density ~ numACGT, data.clustered, group = cluster,
    #cex = log(data.clustered$numACGT / size.threshold, base=10),
    auto.key=T)

xyplot(GCcontent ~ numACGT, data.clustered, group = cluster,
    #cex = log(data.clustered$numACGT / size.threshold, base=10),
    auto.key=T)

Statistics of mitochondrial sequences

which.cluster <- which.min(table(data.clustered$cluster))
data.mt <- subset(data.clustered, cluster == which.cluster, -cluster)
summary(data.mt)

##       name         read.density     GCcontent        numACGT      
##  Min.   : 28457   Min.   : 11.1   Min.   :0.380   Min.   :  2160  
##  1st Qu.:202480   1st Qu.: 24.7   1st Qu.:0.441   1st Qu.:  6998  
##  Median :205008   Median : 27.6   Median :0.448   Median : 19098  
##  Mean   :202409   Mean   : 28.6   Mean   :0.447   Mean   : 26778  
##  3rd Qu.:206727   3rd Qu.: 30.8   3rd Qu.:0.454   3rd Qu.: 34612  
##  Max.   :208054   Max.   :108.7   Max.   :0.495   Max.   :216585

nrow(data.mt)

## [1] 223

sum(data.mt$numACGT)

## [1] 5971387

write.table(data.clustered, 'pg29organelle-scaffolds.fa.clustered.tab',
    quote = FALSE, sep = '\t', row.names = FALSE)