KIRC: comparing GA2x and HiSeq

Plot setup

<<StatisticalPlots`

edfOptions={Frame->True,FrameStyle->Directive[12,(FontFamily->"Helvetica"),Black,Thickness[0.0016]],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,
GridLines->Automatic,GridLinesStyle->Directive[Lighter[Gray]]};

pieChartOptions=BaseStyle->Directive[13,(FontFamily->"Helvetica")]

histogramOptions={Frame->True,FrameStyle->Directive[14,(FontFamily->"Helvetica"),Black,Thickness[0.0016]],
LabelStyle->(FontFamily->"Helvetica"),
GridLines->Automatic,GridLinesStyle->Directive[Lighter[Gray]]};
histogram3DOptions={AxesStyle->Directive[14,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->(FontFamily->"Helvetica")};

listPlotOptions={Frame->True,FrameStyle->Directive[14,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->Directive[14,FontFamily->"Helvetica"],Axes->None,
GridLines->Automatic,GridLinesStyle->Directive[Gray,Opacity[0.8]]};
listPlotOptionsInset={Frame->True,FrameStyle->Directive[20,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,
GridLines->Automatic,GridLinesStyle->Directive[Gray,Opacity[0.6]]};

barchartOptions={Frame->True,BaseStyle->Directive[14,(FontFamily->"Helvetica")],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,FrameTicks->{{Automatic,None},{None,None}},
GridLines->{None,Automatic},GridLinesStyle->Directive[Gray,Opacity[0.4]]};

edfFn[lis_]:=(
y=Table[N[i/Length[lis]],{i,Length[lis]}];
x=Sort[lis];
xy=Table[{x[[i]],y[[i]]},{i,Length[lis]}]
)

(* lis = {{name,value},{name,value},...}*)
edfWithNames[lis_]:=(
sortedLis=Sort[lis,(#1[[2]]<#2[[2]])&];
lisLen=Length[lis];
xy=Table[{sortedLis[[i,1]],sortedLis[[i,2]],N[i/lisLen]},{i,Length[lis]}]
)

edfFn[{1,2,3,4,5}]
edfWithNames[{
{"a",1},
{"b",5},
{"c",3},
{"d",4},
{"e",2}
}
]
BaseStyle->Directive[13,FontFamily->"Helvetica"]
{{1,0.2`},{2,0.4`},{3,0.6`},{4,0.8`},{5,1.`}}
{{"a",1,0.2`},{"e",2,0.4`},{"c",3,0.6`},{"d",4,0.8`},{"b",5,1.`}}

Directories

dataDir="/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples_KIRC/20110531"

KIRC vs KIRP

sampleTypes=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/sample_list.txt","CSV"]

Select[sampleTypes,(StringMatchQ["KIRP",#[[5]]]==True)&]//TableForm

sampleTypeRules=Table[sampleTypes[[i,1]]->sampleTypes[[i,5]],{i,Length[sampleTypes]}]
"MX0281_CGTGAT"/.sampleTypeRules

Driver miRNAs

(*ga2Drivers=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/GaIIdrivers.txt","TSV"];
Dimensions[ga2Drivers]
ga2Drivers[[1;;5,1;;5]]//TableForm*)

HiSeq is not subsampled

(*hiSeqDrivers=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/HiSeqdrivers.txt","TSV"];
Dimensions[hiSeqDrivers]
hiSeqDrivers[[1;;5,1;;5]]//TableForm*)

GA2x data

ga2Tab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples_KIRC/20110531/GaII.raw.txt","TSV"];
ga2Data=Drop[ga2Tab,1];
Dimensions[ga2Data]
ga2Tab[[1;;5,1;;5]]//TableForm

ga2Hdr=ga2Tab[[1]]
ga2HdrA=Table[{i,StringJoin[ga2Hdr[[i]],"_ga"]},{i,Length[ga2Hdr]}]
Dimensions[ga2HdrA]
Dimensions[Transpose[ga2HdrA]]
ga2HdrAT=Transpose[ga2HdrA];

ga2TabA=Drop[Join[ga2HdrAT,ga2Data],1];
Dimensions[ga2TabA]
ga2TabA[[1;;5,1;;5]]//TableForm

hiSeq data

Subsampled data

Clear[hiSeqTab,hiSeqData]
hiSeqTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples_KIRC/20110531/HiSeq.raw.txt","TSV"];
hiSeqData=Drop[hiSeqTab,1];
Dimensions[hiSeqData]
hiSeqTab[[1;;5,1;;5]]//TableForm

Original data, not subsampled

(*hiSeqTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/HiSeq.txt","TSV"];
hiSeqData=Drop[hiSeqTab,1];
Dimensions[hiSeqData]
hiSeqTab[[1;;5,1;;5]]//TableForm*)

hiSeqHdr=hiSeqTab[[1]];
hiSeqHdrA=Table[{i,StringJoin[hiSeqHdr[[i]],"_hs"]},{i,Length[hiSeqHdr]}]

Dimensions[hiSeqHdrA]
Dimensions[Transpose[hiSeqHdrA]]
hiSeqHdrAT=Transpose[hiSeqHdrA];

hiSeqTabA=Drop[Join[hiSeqHdrAT,hiSeqData],1];
Dimensions[hiSeqTabA]
hiSeqTabA[[1;;5,1;;5]]//TableForm

hiSeqTabAData=hiSeqTabA[[All,2;;All]];
{Dimensions[hiSeqTabA],Dimensions[hiSeqTabAData]}
hiSeqTabAData[[1;;5,1;;5]]//TableForm

Compare boxWhisker plots for a particular miRNA

myMiR="hsa-mir-326"

ga2MiRnaRec=Flatten[Select[ga2Data,(StringMatchQ[myMiR,#[[1]]]==True)&]]
hiSeqMiRnaRec=Flatten[Select[hiSeqData,(StringMatchQ[myMiR,#[[1]]]==True)&]]
ga2MiRnaData=Drop[ga2MiRnaRec,1]
hiSeqMiRnaData=Drop[hiSeqMiRnaRec,1]
log10Ga2MiRnadata=Map[Log10,ga2MiRnaData]
log10HiSeqMiRnadata=Map[Log10,hiSeqMiRnaData]

BoxWhiskerChart[{log10Ga2MiRnadata,log10HiSeqMiRnadata},{"Outliers",{"MedianMarker",1,Directive[Thickness[0.006],White]}},ChartLabels->{"GA2x","HiSeq"},BaseStyle->Directive[13,(FontFamily->"Arial")],GridLines->{None,Automatic},GridLinesStyle->Directive[Gray,Opacity[0.5]],FrameLabel->{None,"(Subscript[log, 10])(a)",myMiR},AspectRatio->1.,ImageSize->200]

{Median[hiSeqMiRnaData],Median[ga2MiRnaData],N[Median[hiSeqMiRnaData]/Median[ga2MiRnaData]]}

{Mean[ga2MiRnaData],Variance[ga2MiRnaData]}

joinData=Join[ga2MiRnaData,hiSeqMiRnaData];
{Mean[joinData],Variance[joinData]}

Join the two tables, left-to-right

mainTab=Join[ga2TabA,hiSeqTabAData,2];
Dimensions[mainTab]
mainTab[[1;;5,1;;5]]//TableForm

mainHdr=mainTab[[1]]
mainHdrIds=Drop[mainHdr,1]
mainData=Drop[mainTab,1];
Dimensions[mainData]

Export[ToFileName[dataDir,"expn_matrix.raw.1046-miRs.20110531.txt"],mainTab,"TSV"]

EdgeR labels

labels=Join[ConstantArray[0,56],ConstantArray[1,56]]

forEdgeR=Transpose[{Drop[mainHdr,1],labels}]

Export[ToFileName[dataDir,"labels-for-edgeR.expn_matrix.20110531.csv"],forEdgeR,"CSV"]

Adjusted expression matrix from Elizabeth

(*mainTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/767_nonzero_adjusted/expn_matrix.767-miRs.nonzero-sums.20110516.csv_adjusted.txt","TSV"];
Dimensions[mainTab]
mainTab[[1;;5,1;;5]]//TableForm
mainHdrIds=Drop[mainTab[[1]],1]
*)

Rules

inputIdOrderRules=Table[i->mainHdrIds[[i]],{i,Length[mainHdrIds]}]

Sums

rowTotalFn[row_]:=(Total[row[[2;;All]]])
rowTotalFn[mainTab[[2]]]

zeroRecs=Select[mainData,(rowTotalFn[#]==0)&];
nonZeroRecs=Select[mainData,(rowTotalFn[#]>0)&];
{Length[zeroRecs],Length[nonZeroRecs]}

Records with nonzero row sums

nonZeroTab=Prepend[nonZeroRecs,mainHdr];
Dimensions[nonZeroTab]
nonZeroTab[[1;;5,1;;5]]//TableForm

Export[ToFileName[dataDir,"expn_matrix.768-miRs.nonzero-sums.20110531.txt"],nonZeroTab,"TSV"]

Means

Function

rowMeanFn[row_]:=(Mean[row[[2;;All]]])
rowMeanFn[mainTab[[2]]]

GA2x

ga2Means=Map[rowMeanFn,ga2Data[[1;;All]]];//Timing
ga2Means[[1;;5]]

Ga2 drivers

ga2DriversMeans=Map[rowMeanFn,ga2Drivers[[1;;All]]];//Timing
ga2DriversMeans[[1;;5]]

HiSeq

hiSeqMeans=Map[rowMeanFn,hiSeqData[[1;;All]]];//Timing
hiSeqMeans[[1;;5]]

HiSeq drivers

hiSeqDriversMeans=Map[rowMeanFn,hiSeqDrivers[[1;;All]]];//Timing
hiSeqDriversMeans[[1;;5]]

Combined

means=Map[rowMeanFn,mainTab[[2;;All]]];//Timing
means[[1;;5]]

Variances

edfPlotTitle="KIRC (56): GA2x vs HiSeq"

Function

rowVarianceFn[row_]:=(Variance[row[[2;;All]]])
rowVarianceFn[mainTab[[2]]]

GA2x

ga2Variances=Map[rowVarianceFn,ga2Tab[[2;;All]]];
ga2Variances[[1;;5]]

GA2 drivers

ga2DriversVariances=Map[rowVarianceFn,ga2Drivers[[1;;All]]];
Length[ga2DriversVariances]
ga2DriversVariances[[1;;5]]

HiSeq

hiSeqVariances=Map[rowVarianceFn,hiSeqTab[[2;;All]]];
hiSeqVariances[[1;;5]]

HiSeq drivers

hiSeqDriversVariances=Map[rowVarianceFn,hiSeqDrivers[[1;;All]]];
Length[hiSeqDriversVariances]
hiSeqDriversVariances[[1;;5]]

Combined

variances=Map[rowVarianceFn,mainTab[[2;;All]]];//Timing
variances[[1;;5]]

qVar30=Quantile[variances,0.3]
qVar50=Quantile[variances,0.5]
qVar80=Quantile[variances,0.8]
qVar85=Quantile[variances,0.85]
qVar90=Quantile[variances,0.9]

varPlot=Show[
ListLogLinearPlot[edfFn[variances],Joined->False,PlotRange->{{(10^(-5)),(10^10)},{-0.02,1.02}},PlotStyle->Directive[Darker[Blue],Opacity[0.5],Thickness[0.009]],listPlotOptions,FrameLabel->{"Variance of abundance","Fraction of miRNAs",edfPlotTitle}],

ListLogLinearPlot[edfFn[ga2Variances],Joined->False,PlotStyle->Directive[Darker[Green],Opacity[0.5],Thickness[0.008]]],
ListLogLinearPlot[edfFn[ga2Variances],Joined->False,PlotStyle->Directive[Darker[Red],Opacity[0.5],Thickness[0.007]]],

Graphics[{EdgeForm[{Red,Thickness[0.0015],Opacity[0.5]}],Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],
(*Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar50],1.0}]}],
Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],*)
(*Graphics[{Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],Log[0.003]},{0.0,1.0}]}],*)

ImageSize->350,ImageMargins->5]

aboveQVarianceRecs=Select[mainData,(rowVarianceFn[#]>qVar85)&];
{Dimensions[aboveQVarianceRecs],"miRs"}

{Length[aboveQVarianceRecs[[All,1]]],aboveQVarianceRecs[[All,1]]}

aboveQVarianceRecsPp=Prepend[aboveQVarianceRecs,mainHdr];
aboveQVarianceRecsPp[[1;;4,1;;10]]//TableForm
Dimensions[aboveQVarianceRecsPp]

(*Export[ToFileName[dataDir,"expn_matrix.q85_156-miRs.20110519.txt"],aboveQVarianceRecsPp,"TSV"]*)

Relationship between mean and variance

GA2x

ga2Mv=Transpose[{ga2Means,ga2Variances}];
Dimensions[ga2Mv]
ga2Mv[[1;;5]]//TableForm

ga2DriversMv=Transpose[{ga2DriversMeans,ga2DriversVariances}];
Dimensions[ga2DriversMv]
ga2DriversMv[[1;;5]]//TableForm

ga2MVPlot=Show[
ListLogLogPlot[ga2Mv,Joined->False,PlotRange->{{(10^(-3)),(10^5)},{(10^(-4)),(10^8)}},PlotStyle->Directive[Darker[Green],Opacity[0.8],Thickness[0.009]],listPlotOptions,FrameLabel->{"Variance of abundance","Mean abundance","KIRC: GA2x"}],

ListLogLogPlot[ga2DriversMv,Joined->False,PlotStyle->Directive[Black,Opacity[0.8],Thickness[0.009]]],

(*Graphics[{EdgeForm[{Red,Thickness[0.0015],Opacity[0.5]}],Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar30],1.0}]}],*)
(*Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar50],1.0}]}],
Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],*)
(*Graphics[{Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],Log[0.003]},{0.0,1.0}]}],*)

ImageSize->350,ImageMargins->5,AspectRatio->1]

HiSeq

hiSeqMv=Transpose[{hiSeqMeans,hiSeqVariances}];
Dimensions[hiSeqMv]
hiSeqMv[[1;;5]]//TableForm

hiSeqDriversMv=Transpose[{hiSeqDriversMeans,hiSeqDriversVariances}];
Dimensions[hiSeqDriversMv]
hiSeqDriversMv[[1;;5]]//TableForm

hiSeqMVPlot=Show[
ListLogLogPlot[hiSeqMv,Joined->False,PlotRange->{{(10^(-3)),(10^5)},{(10^(-4)),(10^8)}},PlotStyle->Directive[Darker[Red],Opacity[0.8],Thickness[0.009]],listPlotOptions,FrameLabel->{"Variance of abundance","Mean abundance","KIRC: HiSeq"}],

ListLogLogPlot[hiSeqDriversMv,Joined->False,PlotStyle->Directive[Magenta,Opacity[0.8],Thickness[0.009]]],

(*Graphics[{EdgeForm[{Red,Thickness[0.0015],Opacity[0.5]}],Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar30],1.0}]}],*)
(*Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar50],1.0}]}],
Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],*)
(*Graphics[{Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],Log[0.003]},{0.0,1.0}]}],*)

ImageSize->350,ImageMargins->5,AspectRatio->1]

Combined

{Min[means],Max[means]}

mv=Transpose[{means,variances}];
Dimensions[mv]
mv[[1;;5]]//TableForm

varPlot=Show[
ListLogLogPlot[ga2Mv,Joined->False,PlotRange->{{(10^(-3)),(10^5)},{(10^(-4)),(10^8)}},PlotStyle->Directive[Green,Opacity[0.8],Thickness[0.009]],listPlotOptions,FrameLabel->{"Mean abundance","Variance of abundance",edfPlotTitle}],

ListLogLogPlot[hiSeqMv,Joined->False,PlotStyle->Directive[Red,Opacity[0.6],Thickness[0.009]]],
ListLogLogPlot[ga2DriversMv,Joined->False,PlotStyle->Directive[Black,Opacity[0.8],Thickness[0.009]],PlotMarkers->{△,11}],
ListLogLogPlot[hiSeqDriversMv,Joined->False,PlotStyle->Directive[Black,Opacity[0.8],Thickness[0.009]],PlotMarkers->{◦,11}],

(*Graphics[{EdgeForm[{Red,Thickness[0.0015],Opacity[0.5]}],Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar30],1.0}]}],*)
(*Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar50],1.0}]}],
Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],*)
(*Graphics[{Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],Log[0.003]},{0.0,1.0}]}],*)

Graphics[{Green,Text["GA2x",{Log[0.0015],Log[2*(10^7)]},{-1,0},Background->White]}],
Graphics[{Red,Text["HiSeq",{Log[0.0015],Log[0.5*(10^7)]},{-1,0},Background->White]}],

ImageSize->400,ImageMargins->5,AspectRatio->1,BaseStyle->Directive[12,(FontFamily->"Arial")]]

(*Export["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples/data/differential/GA2x-HiSeq.KIRC.ss1.means-variances.20110525.pdf",varPlot]*)

Colorbar for hclust heatmap

colLabs=ReadList["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/expn_matrix.q30_732-miRs.blue-red.col-labels.20110516.txt","String"];
Length[colLabs]
colLabs[[1;;5]]

colLabsSplit=Table[{
colLabs[[i]],
StringTake[colLabs[[i]],-2]
},{i,Length[colLabs]}];
colLabsSplit[[1;;3]]//TableForm

gaHiSeqColours=Table[
Switch[colLabsSplit[[i,2]],
"ga",Lighter[Green],
"hs",Lighter[Lighter[Red]]
],{i,Length[colLabsSplit]}];

ArrayPlot[{gaHiSeqColours},ImageSize->800,Frame->False,AspectRatio->0.075]

Table[Flatten[{i,colLabsSplit[[i]]}],{i,Length[colLabsSplit]}]//TableForm

Consensus clustering

ccpDir="/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/104_q90/KIRC-56.GA2x-HiSeq.precursor.104-q90.km_spearman_reps500"

clusterOrder={3,5,4,1,2}

png=Import[ToFileName[ccpDir,"consensus.k6.png"],"PNG"]

consensusClass=Import[ToFileName[ccpDir,"KIRC-56.GA2x-HiSeq.precursor.104-q90.km_spearman_reps500.k=6.consensusClass.csv"],"CSV"];
Dimensions[consensusClass]
consensusClass[[1;;5]]//TableForm

ctOrderTab=Import[ToFileName[ccpDir,"cTorder.k6.txt"],"TSV"];
ctOrder=Drop[ctOrderTab,1];
Length[ctOrder]
ctOrder[[1;;5]]//TableForm

ctOrderFull=Table[{
ctOrder[[i,1]],ctOrder[[i,2]],
id=ctOrder[[i,2]]/.inputIdOrderRules,
StringTake[id,13]/.sampleTypeRules,
StringTake[id,-2]
},{i,Length[ctOrder]}];
ctOrderFull//TableForm

gaHiSeqCCPColours=Table[
Switch[ctOrderFull[[i,5]],
"ga",Lighter[Green],
"hs",Lighter[Lighter[Red]]
],{i,Length[ctOrderFull]}];
ArrayPlot[{gaHiSeqCCPColours},ImageSize->800,Frame->False,AspectRatio->0.075]

typeColours=Table[
Switch[ctOrderFull[[i,4]],
"KIRC",Lighter[Orange],
"KIRP",Lighter[Lighter[Blue]]
],{i,Length[ctOrderFull]}];
ArrayPlot[{gaHiSeqCCPColours,typeColours},ImageSize->800,Frame->False,AspectRatio->0.075]

Low-abundance miRNAs

belowQVarianceRecs=Select[mainData,(rowVarianceFn[#]<qVar50)&];
{Dimensions[belowQVarianceRecs],"miRs x columns"}

{Length[belowQVarianceRecs[[All,1]]],"miRNAs>",belowQVarianceRecs[[All,1]]}

belowQVarianceRecsPp=Prepend[belowQVarianceRecs,mainHdr];
belowQVarianceRecsPp[[1;;4,1;;10]]//TableForm
Dimensions[belowQVarianceRecsPp]

nonZeroBelowQVarianceRecs=Select[belowQVarianceRecs,(rowTotalFn[#]>0)&];
{Length[nonZeroBelowQVarianceRecs],Length[belowQVarianceRecs],"miRNAs x columns"}
nonZeroBelowQVarianceRecs[[1;;5,1;;5]]//TableForm
nonZeroBelowQVarianceMiRNAs=nonZeroBelowQVarianceRecs[[All,1]]

mainHdr;
rawSampleIds=Drop[mainHdr,1]
sampleIds=Union[Table[StringTake[rawSampleIds[[i]],12],{i,Length[rawSampleIds]}]]
Length[sampleIds]

nonZeroBelowQVarianceMiRNAs

miR=nonZeroBelowQVarianceMiRNAs[[1]]
sel=Select[mainData,(StringMatchQ[miR,#[[1]]]==True)&];
selPp=Prepend[sel,mainHdr];
selPpT=Drop[Transpose[selPp],1]

baseId=sampleIds[[1]]
selById=Select[selPpT,(StringMatchQ[StringTake[#[[1]],12],baseId]==True)&]
result=Flatten[{baseId,selById[[All,2]]}]
Map[NumberQ,result]

A function for all miRs (i) and all sample basenames (j)

selFn[i_,j_]:=(
miR=nonZeroBelowQVarianceMiRNAs[[i]];
sel=Select[mainData,(StringMatchQ[miR,#[[1]]]==True)&];
selPp=Prepend[sel,mainHdr];
selPpT=Drop[Transpose[selPp],1];

baseId=sampleIds[[j]];
selById=Select[selPpT,(StringMatchQ[StringTake[#[[1]],12],baseId]==True)&];
result=Flatten[{baseId,selById[[All,2]]}]
)
selFn[1,1]
selFn[1,2]
Table[selFn[1,i],{i,Length[sampleIds]}]

t=Flatten[Table[selFn[i,j],{i,Length[nonZeroBelowQVarianceMiRNAs]},{j,Length[sampleIds]}],1];//Timing
t[[1;;5]]

s1=Select[t,(#[[2]]>0&&#[[3]]>0)&];
s1=Sort[s1,(#1[[2]]>#2[[2]])&];
Length[s1]
s1[[1;;25]]//TableForm

ListPlot[Transpose[{s1[[All,2]],s1[[All,3]]}],listPlotOptions,PlotRange->{{0,2},{0,2}},AspectRatio->Automatic,FrameLabel->{"GA2x","HiSeq"}]

Differentially abundant miRNAs using edgeR

edgeRTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples/data/differential/ss5.edgeR.738-miRs.nonzero-sums.20110525.txt","TSV"];
edgeRTab[[1;;5]]//TableForm
edgeRData=Drop[edgeRTab,1];
Dimensions[edgeRData]

logConc=edgeRData[[All,2]];
logFC=edgeRData[[All,3]];

ListPlot[Transpose[{logConc,logFC}],listPlotOptions,PlotRange->{{-30,0},{-3,3}}]

bhPvalThresh=0.05

daMiRNAsEdgeRRecs=Select[edgeRData,(#[[5]]<=bhPvalThresh)&];
Length[daMiRNAsEdgeRRecs]
logConc1=daMiRNAsEdgeRRecs[[All,2]];
logFC1=daMiRNAsEdgeRRecs[[All,3]];
Prepend[daMiRNAsEdgeRRecs,edgeRTab[[1]]]//TableForm

daMiRNAsEdgeR=daMiRNAsEdgeRRecs[[All,1]]

sh=Show[
ListPlot[Transpose[{logConc,logFC}],PlotStyle->Directive[Lighter[Blue],Opacity[0.6],PointSize[0.012]],listPlotOptions,PlotRange->{{-30,0},{-3.5,3.5}},FrameLabel->{"(Subscript[log, 2])(Concentration)","(Subscript[log, 2])(Fold change)"}],
ListPlot[Transpose[{logConc1,logFC1}],PlotStyle->Directive[Red,PointSize[0.014]]],
ImageSize->400,ImageMargins->5
]

Export["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples/data/differential/ss5.edgeR-RCplot.738-miRs.nonzero-sums.20110525.pdf",sh];

2^-20

Table[Flatten[{
recs1[[i]],
2^recs1[[i,2]]
}],{i,Length[recs1]}]//TableForm

DA miRNAs, in edgeR order

daMiRNAMainTabRecs=Select[mainTab,(MemberQ[daMiRNAsEdgeR,#[[1]]]==True)&];
Dimensions[daMiRNAMainTabRecs]

(*daLabs=daMiRNAMainTabRecs[[All,1]]//TableForm
daMeans=Map[rowMeanFn,daMiRNAMainTabRecs];
daVars=Map[rowVarianceFn,daMiRNAMainTabRecs];
daMv=Transpose[{daMeans,daVars}];
SortBy[daMv,First]//TableForm*)

daToPlotTab=Table[
mainRec=First[Select[mainTab,(StringMatchQ[daMiRNAsEdgeR[[i]],#[[1]]]==True)&]];
Flatten[{
daMiRNAsEdgeR[[i]],
rowMeanFn[mainRec],
rowVarianceFn[mainRec]
}],{i,Length[daMiRNAsEdgeR]}];
daToPlotTab//TableForm
daMv=daToPlotTab[[All,2;;3]]

daVarPlot=Show[
ListLogLogPlot[ga2Mv,Joined->False,PlotRange->{{(10^(-3)),(10^5)},{(10^(-4)),(10^8)}},PlotStyle->Directive[Green,Opacity[0.8],Thickness[0.009]],listPlotOptions,FrameLabel->{"Mean abundance","Variance of abundance",edfPlotTitle}],

ListLogLogPlot[hiSeqMv,Joined->False,PlotStyle->Directive[Red,Opacity[0.6],Thickness[0.009]]],
(*ListLogLogPlot[ga2DriversMv,Joined->False,PlotStyle->Directive[Blue,Opacity[0.8],Thickness[0.009]],PlotMarkers->{△,14}],
ListLogLogPlot[hiSeqDriversMv,Joined->False,PlotStyle->Directive[Blue,Opacity[0.8],Thickness[0.009]],PlotMarkers->{◦,14}],*)

ListLogLogPlot[daMv,Joined->False,PlotStyle->Directive[Black,Opacity[1],Thickness[0.009]],PlotMarkers->{□,15}],

Graphics[{EdgeForm[{Blue,Thickness[0.001],Opacity[0.2]}],Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-4))],Log[(10^(-4))]},{Log[0.5],Log[(10^8)]}]}],
(*Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar50],1.0}]}],
Graphics[{EdgeForm[{Red,Thickness[0.002],Opacity[0.5]}],Darker[Blue],Opacity[0.08],Rectangle[{Log[(10^(-5))],0.0},{Log[qVar80],1.0}]}],*)
(*Graphics[{Darker[Blue],Opacity[0.1],Rectangle[{Log[(10^(-5))],Log[0.003]},{0.0,1.0}]}],*)

Graphics[{Green,Text["GA2x",{Log[0.0015],Log[2*(10^7)]},{-1,0},Background->White]}],
Graphics[{Red,Text["HiSeq",{Log[0.0015],Log[0.5*(10^7)]},{-1,0},Background->White]}],

ImageSize->400,ImageMargins->5,AspectRatio->1,BaseStyle->Directive[12,(FontFamily->"Arial")]]

Export["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/HiSeq/subsamples/data/differential//ss5.mean-variance.DA-miRNAs-738.20110525.pdf",daVarPlot]

Write out an expression matrix that has only the DA miRNAs

daMiRNARecsPp=Prepend[daMiRNAMainTabRecs,mainHdr];
Dimensions[daMiRNARecsPp]
daMiRNARecsPp[[All,1;;5]]//TableForm

(*Export[ToFileName[dataDir,"expn_matrix.9-diff-abund-miRs.20110518.txt"],daMiRNARecsPp,"TSV"]*)

Compare DA miRNAs per sample, between GA2x and HiSeq

The two datasets have the same MX IDs

ga2MxIds=Drop[ga2Hdr,1]
hiSeqMxIds=Drop[hiSeqHdr,1]
{Complement[ga2MxIds,hiSeqMxIds],Complement[hiSeqMxIds,ga2MxIds]}

...but the IDs are not in the same order in the two datasets

Table[ga2MxIds[[i]]==hiSeqMxIds[[i]],{i,Length[ga2MxIds]}]

Get data for one miRNA

daMiRNAsEdgeR

myMiRNA="hsa-mir-316"
(*myMiRNA=daMiRNAsEdgeR[[4]]*)

tabSubset=Prepend[
Select[daMiRNARecsPp,(StringMatchQ[myMiRNA,#[[1]]]==True)&],
mainHdr
];
tabSubsetT=Transpose[tabSubset]
tabSubsetTData=Drop[tabSubsetT,1];

getRecsForAMiRFn[mir_]:=(
tabSubset=Prepend[
Select[daMiRNARecsPp,(StringMatchQ[myMiRNA,#[[1]]]==True)&],
mainHdr
];
tabSubsetT=Transpose[tabSubset];
tabSubsetTData=Drop[tabSubsetT,1]
)
myRecs=getRecsForAMiRFn[myMiRNA]

mySample=ga2MxIds[[1]]
Flatten[Select[tabSubsetTData,(StringMatchQ[mySample,StringTake[#[[1]],13]]==True)&]]

getAbundancePairFn[mxId_,recsForAMiR_]:=(
f=Flatten[Select[recsForAMiR,(StringMatchQ[mxId,StringTake[#[[1]],13]]==True)&]];
{mxId,f[[2]],f[[4]]}
)
getAbundancePairFn["MX0281_ACATCG",myRecs]

myMiRNA
myRecs=getRecsForAMiRFn[myMiRNA]

t=Table[getAbundancePairFn[ga2MxIds[[i]],myRecs],{i,Length[ga2MxIds]}];
t=SortBy[t,Last]
Dimensions[t]
Reverse[t]//TableForm

Show[Table[Graphics[{Blue,Arrow[{{t[[i,2]],i},{t[[i,3]],i}}]}],{i,Length[t]}],PlotRange->{{0,12}, All},
Frame->True,FrameLabel->{"Abundance, GA2x->HiSeq","Per sample (ordered by HiSeq)",myMiRNA},FrameStyle->Directive[12,(FontFamily->"Arial"),GridLines->{Automatic,None}],ImageSize->320,AspectRatio->1
]

Parsed result from RNAfold

rnaFoldResult=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/Secondary_structure/tests/mature.hsa.34660-seqs.20-indexes.20110527.fa.rnafold.parsed.txt","TSV"];
(*result=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/batch_effects/Secondary_structure/tests/mature.hsa.constructs_with_CACTGT.fa.rnafold.parsed.txt","TSV"];*)
Dimensions[rnaFoldResult]
rnaFoldResult[[1;;10]]//TableForm

scores=rnaFoldResult[[All,4]];
rsScores=RandomSample[scores,5000];
rsScores[[1;;5]]
freeEnEns=rnaFoldResult[[All,5]];
rsFreeEnEns=RandomSample[freeEnEns,5000];
rsFreeEnEns[[1;;5]]

sortedResult=Sort[rnaFoldResult,(#1[[5]]<#2[[5]])&];
sortedResult[[1;;30]]//TableForm

revSortedResult=Reverse[sortedResult];
revSortedResult[[1;;30]]//TableForm

Abundance on GA2x vs. HiSeq

Abundance data for a miRNA

myMiRNA
myRecs=getRecsForAMiRFn[myMiRNA];

myRecsSplit=Table[
Flatten[{
StringSplit[myRecs[[i,1]],"_"],
myRecs[[i,2;;All]]
}],{i,Length[myRecs]}]

RNAfold results for this miRNA

rnaFoldResult[[1;;3]]//TableForm

s11=Select[rnaFoldResult,(StringMatchQ[myMiRNA,ToLowerCase[#[[1]]]]==True)&];
Length[s11]
s11//TableForm

s11Rules=Table[s11[[i,3]]->s11[[i,4]],{i,Length[s11]}]

ListPlot[edfFn[s11[[All,4]]],listPlotOptions]

myRecsSplit[[1;;5]]//TableForm

myRecsSplitWithFE=Table[
Flatten[{
myRecsSplit[[i]],
myRecsSplit[[i,2]]/.s11Rules
}],{i,Length[myRecsSplit]}];
myRecsSplitWithFE//TableForm

selGA=Select[myRecsSplitWithFE,(StringMatchQ["ga",#[[3]]]==True)&]
selHS=Select[myRecsSplitWithFE,(StringMatchQ["hs",#[[3]]]==True)&]

aGA=selGA[[All,{5,4}]];
aHS=selHS[[All,{5,4}]];

Show[
ListPlot[aGA,PlotRange->{All,{-20,-10}}],
ListPlot[aHS,PlotStyle->Green]
]