Introduction

Let’s QC this data.

library("gplots")
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
library("reshape2")

Multi-qc results

Please have a look at the multiQC report. Here are a few key points:

  • Skewer trimming resulted in loss of only a tiny number of bases. This indicates the sequence quality is very high.

  • Fastqc results showing the number of unique and duplicate reads indicates a few samples with <10M unique reads.

  • Per seqence GC content showed an unusual profile for two samples. PG1423-EOS R1 and R2 had GC profile max at 40% compared to the mean. PG2090-EOS also showed an unusual pattern with underrepresented low GC%.

  • Sequence duplication levels were elevated for some fastq files. Here are the files of concern, with <20% unique reads: PG3627-POD1_S86_R1_001 PG3627-POD1_S86_R2_001 PG3609-T0_S317_R1_001 PG2090-EOS_S134_R1_001 PG2090-EOS_S134_R2_001

  • There were two files with overrepresented sequences: PG2090-EOS R1 and R2. Others are okay.

  • Adapter content was very low which is good.

The fastq files were also checked with validatefastq-assembly which looks for signs of file corruption which can occur in large data transfers. No problematic files were detected.

rRNA amount

Ribosomal RNA carryover can be a source of noise. The proportion should be <10% and there were a few samples in excess of this including PG2020-EOS, PG815-EOS, PG1452-EOS and PG702-POD1.

rrna <- read.table("rrna_stats.txt")
rrna <- rrna[,c(1,5)]
rrna$V1 <- sapply(strsplit(rrna$V1,"\\."),"[[",1)
rrna$V5 <- gsub("\\(","",rrna$V5)
rrna$V5 <- gsub("%","",rrna$V5)
rrna$V5 <- as.numeric(rrna$V5)
str(rrna)
## 'data.frame':    319 obs. of  2 variables:
##  $ V1: chr  "3166-POD1_S266_R1_001" "3166-T0_S265_R1_001" "3167-POD1_S268_R1_001" "3167-T0_S267_R1_001" ...
##  $ V5: num  0.57 1.11 0.61 0.93 0.96 0.79 0.7 5.2 1.14 2.83 ...
rrna2 <- rrna[,2]
names(rrna2) <- rrna[,1]

par(mar=c(5,8,3,1))
barplot(rrna2,horiz=TRUE,las=1,cex.names=0.5,main="rRNA carryover")

rrna2 <- rrna2[order(-rrna2)]
barplot(head(rrna2,20),horiz=TRUE,las=1,cex.names=0.6,main="rRNA carryover")

Load the data

tmp <- read.table("3col.tsv.gz",header=FALSE)
x <- as.matrix(acast(tmp, V2~V1, value.var="V3", fun.aggregate = sum))
x <- as.data.frame(x)
accession <- sapply((strsplit(rownames(x),"\\|")),"[[",2)
symbol<-sapply((strsplit(rownames(x),"\\|")),"[[",6)
x$geneid <- paste(accession,symbol)
xx <- aggregate(. ~ geneid,x,sum)
rownames(xx) <- xx$geneid
colnames <- gsub("T0R","T0",colnames(xx))
xx$geneid = NULL
xx <- round(xx)
xx[1:10,1:6]
##                             3166-POD1 3166-T0 3167-POD1 3167-T0 3171-POD1
## ENSG00000000003.15 TSPAN6           3       1         5       5        23
## ENSG00000000005.6 TNMD              0       0         0       0         0
## ENSG00000000419.14 DPM1           685     577       521     735       811
## ENSG00000000457.14 SCYL3          622     611       550     777       789
## ENSG00000000460.17 C1orf112       181     171       232     263       215
## ENSG00000000938.13 FGR          33797   44344     31524   38959     26402
## ENSG00000000971.16 CFH            106      40        98     183       195
## ENSG00000001036.14 FUCA2         1229     769      1150     868       978
## ENSG00000001084.13 GCLC           944    1085       577     961       908
## ENSG00000001167.15 NFYA          1243    1277      1295    1605      1166
##                             3171-T0
## ENSG00000000003.15 TSPAN6         4
## ENSG00000000005.6 TNMD            1
## ENSG00000000419.14 DPM1         494
## ENSG00000000457.14 SCYL3        575
## ENSG00000000460.17 C1orf112     196
## ENSG00000000938.13 FGR        33751
## ENSG00000000971.16 CFH          130
## ENSG00000001036.14 FUCA2        805
## ENSG00000001084.13 GCLC         798
## ENSG00000001167.15 NFYA        1251

Number of reads per sample

Let’s look at the number of reads per sample

Most samples were in the range of 25-30 million assigned reads. Just 2 samples had less than 20 million reads: PG1452-EOS and PG1423-EOS. The maximum read count was about 40 million for PG7072-EOS.

xxcs <- colSums(xx)
par(mar=c(5,8,3,1))
barplot(xxcs,horiz=TRUE,las=1,main="no. reads per sample")