<p>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 <a href="http://bmcmolbiol.biomedcentral.com/articles/10.1186/1471-2199-7-3" rel="nofollow noreferrer">RIN numbers</a> and <a href="http://www.agilent.com/cs/library/applications/5991-5258EN.pdf" rel="nofollow noreferrer">DIN numbers</a>, with a description of the variables that are included in the model.</p>

<p>The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.</p>

<p>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:</p>

<pre><code>myfiles&lt;-list.files()
results&lt;-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n &lt;- which(substr(myfiles,69,79)=="Results.csv")
targets &lt;- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2&lt;- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3&lt;-paste0(substr(targets2[1,2],1,67),"_S")
chip_num&lt;-length(n)
samp_num_chip&lt;-length(grep(str3,myfiles))
samp_num&lt;-length(which(substr(myfiles,69,74)=="Sample"))
samples1&lt;- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
    ##mod chip numb
    targets2&lt;- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
    str3&lt;-paste0(substr(targets2[1,2],1,67),"_S")
    samp_num_chip&lt;-length(grep(str3,myfiles))
    ##mod sample numb
    str1&lt;-substr(targets2[1,2],1,67)
    for(j in 1:samp_num_chip)
    {
        iii&lt;-j
        str5&lt;-paste0(str1,"_Sample",iii,".csv")
        samples1[j,i]&lt;-str5
    }
}
samples1&lt;-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 &lt;- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat &lt;- matrix(NA, ncol=1, nrow=nrow(elec))
dta &lt;- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
    x &lt;- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
    dta[,i] &lt;- x[,2]
}
time&lt;- x[,1]
for(samples in 1:length(samples1))
{
    ##store max fluor in max.peaks matrix
    lad&lt;-samples
    ti1&lt;-17
    ti2&lt;-65
    p.row&lt;-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
    peaks &lt;- matrix(NA, nrow=p.row, ncol=1)
    max.peaks &lt;- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
    for (i in 1:length(seq(ti1,ti2,.05)))
    {
        ii&lt;-seq(ti1,ti2,.05)[i]
        x0&lt;- ii
        x0&lt;-round(x0,digits=3)
        y0&lt;- ii
        x&lt;- ii+.05
        x&lt;-round(x,digits=3)
        y&lt;- ii
        max.peaks[i]&lt;- max(dta[which(time==x0):which(time==x),lad])
    }    
    ##human or mouse
    cent&lt;-max(max.peaks)*.6
    ##flatworm
    ##cent&lt;-max(max.peaks)*.05
    ##define the middle of the segments
    max.avg&lt;-seq(ti1+.025,ti2+.025,.05)
    ##find time associated with peaks &gt; 5 fu and fill seconds matrix with them
    ##12.5% of max fu incase chip is fu
    sizes &lt;- which(max.peaks&gt;cent)
    gt35&lt;-which(max.avg[sizes]&gt;35)
    lt55&lt;-which(max.avg[sizes]&lt;55)
    tfr&lt;-which(duplicated(c(gt35,lt55))==TRUE)
    area&lt;-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&lt;-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&lt;-time2[which(area[,1]&lt;cent),1][1]
    t2&lt;-time2[which(area[,1]&lt;cent),1][2]
    t3&lt;-time2[which(area[,3]&lt;cent),3][1]
    t4&lt;-time2[which(area[,3]&lt;cent),3][2]
    s18a&lt;-ifelse(is.na(t1,-30,round(t1,digits=1)))
    s18b&lt;-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
    s28a&lt;-ifelse(is.na(t2,-30,round(t2,digits=1)))
    s28b&lt;-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
    eta&lt;-which(time==s18a)
    etb&lt;-which(time==s18b)
    tea&lt;-which(time==s28a)
    teb&lt;-which(time==s28b)
    s18&lt;-sum(dta[eta:etb,lad])
    s28&lt;-sum(dta[tea:teb,lad])
    qc.mat[samples]&lt;-  (-1*exp((s28/s18)*-1)+1)*10
}
</code></pre>

<p>Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.</p>

CC BY-SA 3.0
2017-07-17T02:10:55.023Z
2017-07-17T02:10:55.023Z
gringer
<p>The best I could find by looking at github.com is:</p>

<p><a href="https://github.com/brianfleharty/Agilent-BioAnalyzer-RIN-R-code-for-Fly/blob/master/RNA_RIN.R" rel="nofollow noreferrer">https://github.com/brianfleharty/Agilent-BioAnalyzer-RIN-R-code-for-Fly/blob/master/RNA_RIN.R</a></p>

<hr>

<pre class="lang-r prettyprint-override"><code>#set the working directory
dir&lt;- "H:\\B_RIN\\emd"
setwd(dir)
myfiles&lt;-list.files()
which(substr(myfiles,69,79)=="Results.csv")
results&lt;-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n &lt;- which(substr(myfiles,69,79)=="Results.csv")
targets &lt;- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
length(n)
targets2&lt;- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str1&lt;-substr(targets2[1,2],1,67)
str1
str2&lt;-"_S"
str3&lt;-paste(str1,str2,sep="")
str3
chip_num&lt;-length(n)
samp_num_chip&lt;-length(grep(str3,myfiles))
samp_num&lt;-length(which(substr(myfiles,69,74)=="Sample"))
samples1&lt;- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
#mod chip numb
targets2&lt;- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
str1&lt;-substr(targets2[1,2],1,67)
str2&lt;-"_S"
str3&lt;-paste(str1,str2,sep="")
samp_num_chip&lt;-length(grep(str3,myfiles))
#mod sample numb
str1&lt;-substr(targets2[1,2],1,67)
str4&lt;-"_Sample"
for(j in 1:samp_num_chip)
{
iii&lt;-j
end.st&lt;-".csv"
str5&lt;-paste(str1,str4,iii,end.st,sep="")
samples1[j,i]&lt;-str5
}
}
samples1
grep("2100",samples1)
samples1&lt;-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 &lt;- colorRamp(c("darkmagenta","white"))
elec &lt;- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat &lt;- matrix(NA, ncol=1, nrow=nrow(elec))
dta &lt;- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
x &lt;- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
dta[,i] &lt;- x[,2]
}
time&lt;- 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] &gt; 5)
#draw some segments and store max fluor in max.peaks matrix
lad&lt;-samples
ti1&lt;-17
ti2&lt;-65
p.row&lt;-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
peaks &lt;- matrix(NA, nrow=p.row, ncol=1)
max.peaks &lt;- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
for (i in 1:length(seq(ti1,ti2,.05)))
{
ii&lt;-seq(ti1,ti2,.05)[i]
x0&lt;- ii
x0&lt;-round(x0,digits=3)
y0&lt;- ii
x&lt;- ii+.05
x&lt;-round(x,digits=3)
y&lt;- ii
#segments(x0,y0,x,y)
max.peaks[i]&lt;- max(dta[which(time==x0):which(time==x),lad])
}    
#human or mouse
cent&lt;-max(max.peaks)*.6
#flatworm
#cent&lt;-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&lt;-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 &gt; 5 fu and fill seconds matrix with them
#12.5% of max fu incase chip is fu
p1&lt;-which(max.avg==38)
p2&lt;-which(max.avg==55)
sizes &lt;- which(max.peaks&gt;cent)
max.avg[sizes]
gt35&lt;-which(max.avg[sizes]&gt;35)
lt55&lt;-which(max.avg[sizes]&lt;55)
#as.numeric(duplicated(c(gt35,lt55)))
tfr&lt;-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&lt;-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&lt;-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&lt;-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&lt;-area[which(area[,1]&lt;cent),1][1]
t1&lt;-time2[which(area[,1]&lt;cent),1][1]
p2&lt;-area[which(area[,1]&lt;cent),1][2]
t2&lt;-time2[which(area[,1]&lt;cent),1][2]
#area[which(area[,2]&lt;cent),2]
#area[which(area[,2]&lt;cent),2]
p3&lt;-area[which(area[,3]&lt;cent),3][1]
t3&lt;-time2[which(area[,3]&lt;cent),3][1]
p4&lt;-area[which(area[,3]&lt;cent),3][2]
t4&lt;-time2[which(area[,3]&lt;cent),3][2]
#abline(v=area[which(area[,3]&lt;cent),3])
#abline(h=11)
#abline(v=42.75)
v1&lt;-t1
v2&lt;-t2
v3&lt;-t3
v4&lt;-t4
abline(v=v1)
abline(v=v2)
abline(v=v3)
abline(v=v4)
#time[area[which(area[,1]&lt;cent),1][1]==dta[,lad]]
#time[area[which(area[,3]&lt;cent),3][1]==dta[,lad]]
#time[area[which(area[,1]&lt;cent),1][2]==dta[,lad]]
#time[area[which(area[,3]&lt;cent),3][2]==dta[,lad]]
#abline(v=time[area[which(area[,1]&lt;cent),1][1]==dta[,lad]])
#abline(v=time[area[which(area[,3]&lt;cent),3][1]==dta[,lad]])
#abline(v=time[area[which(area[,1]&lt;cent),1][2]==dta[,lad]])
#abline(v=time[area[which(area[,3]&lt;cent),3][2]==dta[,lad]])
#abline(h=area[which(area[,3]&lt;cent),3][1])
#abline(h=area[which(area[,3]&lt;cent),3][2]+3)
s18a&lt;-round(t1,digits=1)
s18b&lt;-round(t3,digits=1)
s28a&lt;-round(t2,digits=1)
s28b&lt;-round(t4,digits=1)
if(is.na(s18a)){s18a&lt;-30}
if(is.na(s18b)){s18b&lt;-35.05}
if(is.na(s28a)){s28a&lt;-30}
if(is.na(s28b)){s28b&lt;-30.05}
eta&lt;-which(time==s18a)
etb&lt;-which(time==s18b)
tea&lt;-which(time==s28a)
teb&lt;-which(time==s28b)
#dta[eta:etb,lad]
#dta[tea:teb,lad]
s18&lt;-sum(dta[eta:etb,lad])
s28&lt;-sum(dta[tea:teb,lad])
s28/s18
qc.mat[samples]&lt;-  (-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 &lt;- read.delim("results-files.txt", sep=",",header=F,skip=0,as.is=T)
#targets &lt;- read.delim(results[1,], sep=",",header=F,skip=14,as.is=T)
#targets[1:5,]
tbl&lt;-which(targets[,1]=="Overall Results:")
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples&lt;-as.character(targets[which(targets[,1]=="Sample Name"),2])
#read in the targets file                               
colnms&lt;-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
mat1 &lt;- matrix(NA, ncol=length(colnms), nrow=length(samples)*length(results))
colnames(mat1)&lt;-colnms
rownames(mat1)&lt;-rep(samples,length(results))
for(i in 1:length(results)) 
{
targets &lt;- 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&lt;-as.character(targets[which(targets[,1]=="Sample Name"),2])
tbl&lt;-which(targets[,1]=="Overall Results:")
length(which(targets[,1]=="Overall Results:"))
tbl[1]
tbl[1]
colnms&lt;-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
for (j in 1:length(samples))
{
mat1[((length(samples)*(i-1))+j),]&lt;-as.matrix(targets[(tbl[j]+1):(tbl[j]+4),2])
}
}
mat1
cbind(mat1,qc.mat)
colnames(qc.mat)&lt;-"B-RIN"
write.table(cbind(mat1,qc.mat),file=paste(dir,"\\RIN_VS_B-RIN.txt",sep=""), sep="\t",col.names=NA)
</code></pre>

CC BY-SA 3.0
2017-07-12T10:23:40.667Z
Thanks for posting and answering your own question! It would be great if you could also explain how to use this code to do what your question asks for.
terdon
CC BY-SA 3.0
2017-07-12T09:38:45.987Z
2017-07-12T10:23:21.807Z
2017-07-12T10:23:21.807Z
719016
<p>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.</p>

<p>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?</p>

<phylogenetics><quantification>
how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation?
# how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation?

(Asked by: username_0 on 2017-07-12)

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?

**Tags:** phylogenetics, quantification

---

## 2 answers

### Answer by username_1 on 2017-07-17. Score: 2

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]
Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.


---

### Answer by username_0 on 2017-07-12. Score: 1

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]

> **Comments:**
> **username_2**: Thanks for posting and answering your own question! It would be great if you could also explain how to use this code to do what your question asks for.  

---

