1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
|
package align2;
import java.util.ArrayList;
import java.util.Locale;
import fileIO.ReadWrite;
import stream.ConcurrentLegacyReadInputStream;
import stream.RTextInputStream;
import stream.Read;
import stream.SiteScore;
import structures.ListNum;
public class MakeQualityHistogram {
public static void main(String[] args){
String fname1=args[0];
String fname2=(args.length>1 ? args[1] : null);
assert(fname2==null || !fname1.equals(fname2)) : "Error - input files have same name.";
long maxReads=0;
RTextInputStream rtis=new RTextInputStream(fname1, fname2, maxReads);
ConcurrentLegacyReadInputStream cris=new ConcurrentLegacyReadInputStream(rtis, maxReads);
int[][][] counts=process(cris);
printMappedHistogram(counts[0]);
System.out.println();
printPairedHistogram(counts[1]);
// System.out.println("*** main() finished ***");
}
public static void printMappedHistogram(int[][] mapped){
System.out.println("#Error Quality Histogram");
System.out.println("Quality\tMapped\tUnmapped\tPercent Mapped");
for(int i=0; i<mapped[0].length; i++){
int e=mapped[0][i];
int m=mapped[1][i];
float percent=e*100f/(e+m);
System.out.println(i+"\t"+e+"\t"+m+"\t"+String.format(Locale.ROOT, "%.3f", percent));
}
}
public static void printPairedHistogram(int[][] paired){
System.out.println("#Error Quality Histogram");
System.out.println("Quality\tPaired\tSingle\tPercent Paired");
for(int i=0; i<paired[0].length; i++){
int e=paired[0][i];
int m=paired[1][i];
float percent=e*100f/(e+m);
System.out.println(i+"\t"+e+"\t"+m+"\t"+String.format(Locale.ROOT, "%.3f", percent));
}
}
public static int[][][] process(ConcurrentLegacyReadInputStream cris){
cris.start();
int[][] mapped=new int[2][50];
int[][] paired=new int[2][50];
ListNum<Read> ln=cris.nextList();
ArrayList<Read> readlist=ln.list;
while(!readlist.isEmpty()){
processList(readlist, mapped, paired);
cris.returnList(ln.id, readlist.isEmpty());
//System.err.println("Waiting on a list...");
ln=cris.nextList();
readlist=ln.list;
}
//System.err.println("Returning a list... (final)");
assert(readlist.isEmpty());
cris.returnList(ln.id, readlist.isEmpty());
ReadWrite.closeStream(cris);
return new int[][][] {mapped, paired};
}
private static void processList(ArrayList<Read> list, int[][] mapped, int[][] paired) {
for(Read r : list){
processRead(r, mapped, paired);
// if(r.mate!=null){
// processRead(r.mate, mapped, paired);
// }
}
}
private static void processRead(Read r, int[][] mapped, int[][] paired) {
if(r.chrom<1 && r.numSites()>0){
SiteScore ss=r.topSite(); //Should not be necessary
r.start=ss.start;
r.stop=ss.stop;
r.chrom=ss.chrom;
r.setStrand(ss.strand);
}
int avgQ=r.avgQualityInt(true, 0);
if(r.chrom>0){
mapped[0][avgQ]++;
}else{
mapped[1][avgQ]++;
}
if(r.paired()){
paired[0][avgQ]++;
}else{
paired[1][avgQ]++;
}
}
}
|
