For DNA/RNA quantification machines like the Bioanalyzer or TapeStation, the DNA Integrity Number (DIN) or RNA Integrity Number (RIN) numbers are quoted as a measure of the fragmentation of the material.
How is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation? Is there any open-source software to parse/collate the data from the machines?
Note that DIN (DNA Integrity Number) and RIN (RNA Integrity Number) are different scores. I've had trouble in the past finding code or formulas for calculating both of these scores, which is very frustrating considering how frequently they are used in research papers and for NGS QC. The closest that I have been able to get to this are demonstrative graphs for various RIN numbers and DIN numbers, with a description of the variables that are included in the model.
The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.
Just in case it's helpful, here is an extract of that code, only including lines that seem to be directly involved in calculating the actual B_RIN values. This is not runnable, because it depends on particular input files, and it's still pretty unintelligible:
myfiles<-list.files()
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
##mod chip numb
targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
samp_num_chip<-length(grep(str3,myfiles))
##mod sample numb
str1<-substr(targets2[1,2],1,67)
for(j in 1:samp_num_chip)
{
iii<-j
str5<-paste0(str1,"_Sample",iii,".csv")
samples1[j,i]<-str5
}
}
samples1<-samples1[grep("2100",samples1)]
##RIN Calculator#######refined peakfinder####################
##read bioanalyzer data into a matrix called dta
##since the total RNA and mRNA assays run differently
##skip more lines for the total RNA assay
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
dta[,i] <- x[,2]
}
time<- x[,1]
for(samples in 1:length(samples1))
{
##store max fluor in max.peaks matrix
lad<-samples
ti1<-17
ti2<-65
p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
peaks <- matrix(NA, nrow=p.row, ncol=1)
max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
for (i in 1:length(seq(ti1,ti2,.05)))
{
ii<-seq(ti1,ti2,.05)[i]
x0<- ii
x0<-round(x0,digits=3)
y0<- ii
x<- ii+.05
x<-round(x,digits=3)
y<- ii
max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
}
##human or mouse
cent<-max(max.peaks)*.6
##flatworm
##cent<-max(max.peaks)*.05
##define the middle of the segments
max.avg<-seq(ti1+.025,ti2+.025,.05)
##find time associated with peaks > 5 fu and fill seconds matrix with them
##12.5% of max fu incase chip is fu
sizes <- which(max.peaks>cent)
gt35<-which(max.avg[sizes]>35)
lt55<-which(max.avg[sizes]<55)
tfr<-which(duplicated(c(gt35,lt55))==TRUE)
area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],
max.peaks[sizes][c(gt35,lt55)[tfr]],
max.peaks[sizes+1][c(gt35,lt55)[tfr]])
time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],
max.avg[sizes][c(gt35,lt55)[tfr]],
max.avg[sizes+1][c(gt35,lt55)[tfr]])
t1<-time2[which(area[,1]<cent),1][1]
t2<-time2[which(area[,1]<cent),1][2]
t3<-time2[which(area[,3]<cent),3][1]
t4<-time2[which(area[,3]<cent),3][2]
s18a<-ifelse(is.na(t1,-30,round(t1,digits=1)))
s18b<-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
s28a<-ifelse(is.na(t2,-30,round(t2,digits=1)))
s28b<-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
eta<-which(time==s18a)
etb<-which(time==s18b)
tea<-which(time==s28a)
teb<-which(time==s28b)
s18<-sum(dta[eta:etb,lad])
s28<-sum(dta[tea:teb,lad])
qc.mat[samples]<- (-1*exp((s28/s18)*-1)+1)*10
}
Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.
The best I could find by looking at github.com is:
https://github.com/brianfleharty/Agilent-BioAnalyzer-RIN-R-code-for-Fly/blob/master/RNA_RIN.R
#set the working directory
dir<- "H:\\B_RIN\\emd"
setwd(dir)
myfiles<-list.files()
which(substr(myfiles,69,79)=="Results.csv")
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
length(n)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str1
str2<-"_S"
str3<-paste(str1,str2,sep="")
str3
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
#mod chip numb
targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str2<-"_S"
str3<-paste(str1,str2,sep="")
samp_num_chip<-length(grep(str3,myfiles))
#mod sample numb
str1<-substr(targets2[1,2],1,67)
str4<-"_Sample"
for(j in 1:samp_num_chip)
{
iii<-j
end.st<-".csv"
str5<-paste(str1,str4,iii,end.st,sep="")
samples1[j,i]<-str5
}
}
samples1
grep("2100",samples1)
samples1<-samples1[grep("2100",samples1)]
samples1
write.table(samples1,file=paste(dir,"\\ba_lane.txt",sep=""), sep="\t",col.names=NA)
#RIN Calculator#######refined peakfinder####################
#read bioanalyzer data into a matrix called dta
#since the total RNA and mRNA assays run differently
#skip more lines for the total RNA assay
ramp <- colorRamp(c("darkmagenta","white"))
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
dta[,i] <- x[,2]
}
time<- x[,1]
#plot ladder and define peak threshold in HD
par(mfrow=c(3,4))
for(samples in 1:length(samples1))
#for(samples in 1:nrow(elec))
{
#which(dta[,lad] > 5)
#draw some segments and store max fluor in max.peaks matrix
lad<-samples
ti1<-17
ti2<-65
p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
peaks <- matrix(NA, nrow=p.row, ncol=1)
max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
for (i in 1:length(seq(ti1,ti2,.05)))
{
ii<-seq(ti1,ti2,.05)[i]
x0<- ii
x0<-round(x0,digits=3)
y0<- ii
x<- ii+.05
x<-round(x,digits=3)
y<- ii
#segments(x0,y0,x,y)
max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
}
#human or mouse
cent<-max(max.peaks)*.6
#flatworm
#cent<-max(max.peaks)*.05
plot(time,dta[,lad],type="h",col="dodgerblue",xlab="seconds",ylab="Fluorescence Units",xlim=c(17,75),ylim=c(0,(max(max.peaks)*1.2)),main=elec[lad,1])
abline(h=cent,col="grey")
#define the middle of the segments
max.avg<-seq(ti1+.025,ti2+.025,.05)
#plot the maximum value of each segment
max.peaks
lines(max.avg,max.peaks,type="b",col="deepskyblue")
#lines(max.avg,min.peaks,type="b",col="blue")
#find time associated with peaks > 5 fu and fill seconds matrix with them
#12.5% of max fu incase chip is fu
p1<-which(max.avg==38)
p2<-which(max.avg==55)
sizes <- which(max.peaks>cent)
max.avg[sizes]
gt35<-which(max.avg[sizes]>35)
lt55<-which(max.avg[sizes]<55)
#as.numeric(duplicated(c(gt35,lt55)))
tfr<-which(duplicated(c(gt35,lt55))==TRUE)
c(gt35,lt55)[tfr]
max.avg[sizes][c(gt35,lt55)[tfr]]
max.peaks[sizes][c(gt35,lt55)[tfr]]
points(max.avg[sizes+1][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes-1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],col="deeppink3",pch=19)
#colnms2<-c("max.avg[sizes]","max.peaks[sizes]","max.avg[sizes+1]","max.peaks[sizes+1]","max.avg[sizes-1]")
max.avg[sizes]
max.peaks[sizes]
max.avg[sizes+1]
max.peaks[sizes+1]
max.avg[sizes-1]
max.peaks[sizes-1]
area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]])
time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.avg[sizes][c(gt35,lt55)[tfr]],max.avg[sizes+1][c(gt35,lt55)[tfr]])
p1<-area[which(area[,1]<cent),1][1]
t1<-time2[which(area[,1]<cent),1][1]
p2<-area[which(area[,1]<cent),1][2]
t2<-time2[which(area[,1]<cent),1][2]
#area[which(area[,2]<cent),2]
#area[which(area[,2]<cent),2]
p3<-area[which(area[,3]<cent),3][1]
t3<-time2[which(area[,3]<cent),3][1]
p4<-area[which(area[,3]<cent),3][2]
t4<-time2[which(area[,3]<cent),3][2]
#abline(v=area[which(area[,3]<cent),3])
#abline(h=11)
#abline(v=42.75)
v1<-t1
v2<-t2
v3<-t3
v4<-t4
abline(v=v1)
abline(v=v2)
abline(v=v3)
abline(v=v4)
#time[area[which(area[,1]<cent),1][1]==dta[,lad]]
#time[area[which(area[,3]<cent),3][1]==dta[,lad]]
#time[area[which(area[,1]<cent),1][2]==dta[,lad]]
#time[area[which(area[,3]<cent),3][2]==dta[,lad]]
#abline(v=time[area[which(area[,1]<cent),1][1]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][1]==dta[,lad]])
#abline(v=time[area[which(area[,1]<cent),1][2]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][2]==dta[,lad]])
#abline(h=area[which(area[,3]<cent),3][1])
#abline(h=area[which(area[,3]<cent),3][2]+3)
s18a<-round(t1,digits=1)
s18b<-round(t3,digits=1)
s28a<-round(t2,digits=1)
s28b<-round(t4,digits=1)
if(is.na(s18a)){s18a<-30}
if(is.na(s18b)){s18b<-35.05}
if(is.na(s28a)){s28a<-30}
if(is.na(s28b)){s28b<-30.05}
eta<-which(time==s18a)
etb<-which(time==s18b)
tea<-which(time==s28a)
teb<-which(time==s28b)
#dta[eta:etb,lad]
#dta[tea:teb,lad]
s18<-sum(dta[eta:etb,lad])
s28<-sum(dta[tea:teb,lad])
s28/s18
qc.mat[samples]<- (-1*exp((s28/s18)*-1)+1)*10
text(60,cent*.75,c(format(qc.mat[samples],digits=2),"\n\nB-RIN"))
}
qc.mat
###end of ratio calc
###
###
#barplot(qc.mat,beside=T,ylim=c(0,11))
#set the working directory
#results <- read.delim("results-files.txt", sep=",",header=F,skip=0,as.is=T)
#targets <- read.delim(results[1,], sep=",",header=F,skip=14,as.is=T)
#targets[1:5,]
tbl<-which(targets[,1]=="Overall Results:")
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
#read in the targets file
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
mat1 <- matrix(NA, ncol=length(colnms), nrow=length(samples)*length(results))
colnames(mat1)<-colnms
rownames(mat1)<-rep(samples,length(results))
for(i in 1:length(results))
{
targets <- read.delim(results[i], sep=",",header=F,skip=14,as.is=T)
targets[1:5,]
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
tbl<-which(targets[,1]=="Overall Results:")
length(which(targets[,1]=="Overall Results:"))
tbl[1]
tbl[1]
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
for (j in 1:length(samples))
{
mat1[((length(samples)*(i-1))+j),]<-as.matrix(targets[(tbl[j]+1):(tbl[j]+4),2])
}
}
mat1
cbind(mat1,qc.mat)
colnames(qc.mat)<-"B-RIN"
write.table(cbind(mat1,qc.mat),file=paste(dir,"\\RIN_VS_B-RIN.txt",sep=""), sep="\t",col.names=NA)
