Analyze SNVs

Set up plots

<<StatisticalPlots`

histogramOptions={Frame->True,FrameStyle->Directive[12,(FontFamily->"Helvetica"),Black,Thickness[0.0016]],
LabelStyle->(FontFamily->"Helvetica"),
GridLines->Automatic,GridLinesStyle->Directive[Lighter[Gray]]};
edfOptions={Frame->True,FrameStyle->Directive[12,(FontFamily->"Helvetica"),Black,Thickness[0.0016]],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,
GridLines->Automatic,GridLinesStyle->Directive[Lighter[Gray]]};
histogram3DOptions={AxesStyle->Directive[12,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->(FontFamily->"Helvetica")};
listPlotOptions={Frame->True,FrameStyle->Directive[12,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,
GridLines->Automatic,GridLinesStyle->Directive[Gray,Opacity[0.6]]};
barchartOptions={Frame->True,FrameStyle->Directive[10,(FontFamily->"Helvetica"),Black,Thickness[0.001]],
LabelStyle->(FontFamily->"Helvetica"),Axes->None,FrameTicks->{{Automatic,None},{Automatic,None}},
GridLines->{None,Automatic},GridLinesStyle->Directive[Gray,Opacity[0.4]]};
barchartOptionsA={BaseStyle->Directive[11,FontFamily->"Arial"],GridLines->{Automatic,None},GridLinesStyle->Directive[Gray,Opacity[0.4]]};

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

WashU genomic patient vs library IDs

Normal genomic samples

normalWashUTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/AML/data/genomic_wash-u/washu_AML_N.txt","TSV"];
normalWashUData=Drop[normalWashUTab,1];
Dimensions[normalWashUData]
normalWashUData[[1;;3]]

normalWashULib2PatientRules=Table[normalWashUData[[i,1]]->normalWashUData[[i,2]],{i,Length[normalWashUData]}]
normalWashULibIDs=normalWashUData[[All,1]]
normalWashUPatientIDs=normalWashUData[[All,2]]

Tumor genomic samples

tumorWashUTab=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/AML/data/genomic_wash-u/washu_AML_T.txt","TSV"];
tumorWashUData=Drop[tumorWashUTab,1];
Dimensions[tumorWashUData]
tumorWashUData[[1;;3]]

tumorWashULib2PatientRules=Table[tumorWashUData[[i,2]]->tumorWashUData[[i,1]],{i,Length[tumorWashUData]}]
tumorWashULibIDs=tumorWashUData[[All,1]]
tumorWashUPatientIDs=tumorWashUData[[All,2]]

Check that all pairs are paired

Intersection[normalWashUPatientIDs,tumorWashUPatientIDs]
Complement[normalWashUPatientIDs,tumorWashUPatientIDs]

Combine the tables

washUdata=Join[normalWashUData,tumorWashUData];
washUdata2col=washUdata[[All,1;;2]]
{Length[Union[washUdata2col[[All,2]]]],Union[washUdata2col[[All,2]]]}

Patient IDs from Wash U

patientIdsWashU=Union[washUdata2col[[All,2]]]

AML subtype information

Import TSV fle

subtypeInfo=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/AML/data/subtypes/051110 UPN linked to TCGA number with cyto.txt","TSV"];
subtypeHdr=subtypeInfo[[1]]
subtypeData=Drop[subtypeInfo,1];
Dimensions[subtypeData]
subtypeInfo[[1;;4]]//TableForm

Get all patient IDs

patientIdsSubtypeData=subtypeData[[All,3]]
uniqPatientIdsSubtypeData=Union[patientIdsSubtypeData]
{Length[patientIdsSubtypeData],Length[uniqPatientIdsSubtypeData]}

Compare patient IDs from subtype data and WashU genomic data

{Length[Intersection[patientIdsSubtypeData,patientIdsWashU]],Intersection[patientIdsSubtypeData,patientIdsWashU]}
Complement[patientIdsWashU,patientIdsSubtypeData]
Complement[patientIdsSubtypeData,patientIdsWashU]

Get all library IDs

rnaseqLibIds=subtypeData[[2;;All,1]]
Length[rnaseqLibIds]

Get and count subtypes

mOrder={"M0","M1","M2","M3","M4","M5","M6","M7","nc","?"}

fabIds=subtypeData[[All,9]] ;
tallyFAB=Tally[fabIds];
sTallyFAB=Sort[tallyFAB,(#1[[2]]>#2[[2]])&];
sTallyFAB//TableForm
BarChart[Reverse[sTallyFAB[[All,2]]],BarOrigin->Left,ChartLabels->Reverse[sTallyFAB[[All,1]]],barchartOptionsA,ImageSize->350,AspectRatio->1/GoldenRatio,AxesLabel->{"Samples","Subtype"}]

subtypeLibraryIds=subtypeInfo[[2;;All,{1,9}]]

Rules that link RNAseq library ID and patient ID

rnaSeqLibId2PatientIDRules=Table[subtypeData[[i,1]]->subtypeData[[i,3]],{i,Length[subtypeData]}];
patientID2rnaSeqLibIdRules=Table[subtypeData[[i,3]]->subtypeData[[i,1]],{i,Length[subtypeData]}];
{Length[rnaSeqLibId2PatientIDRules],"rules"}
"A00002"/.rnaSeqLibId2PatientIDRules
2807/.patientID2rnaSeqLibIdRules

All rules

SNV data from the VarDB

Import

snvs=Import["/Users/grobertson/GENEREG/ChIP/CHIP_SEQ/SOLEXA/TCGA/AML/washu_genomic/AML/AML_SNVs.txt","TSV"];
Dimensions[snvs]

snvs[[1;;5]]//TableForm

Count how many times each SNV (coord) occurs

tallyIds=Tally[snvs[[All,2]]];
Length[tallyIds]
sortedTallyIds=Sort[tallyIds,(#1[[2]]>#2[[2]])&];
tallyPlot=ListLogLinearPlot[edf[sortedTallyIds[[All,2]]],PlotRange->{All,All},listPlotOptions,Joined->True,PlotStyle->Directive[Thickness[0.008]],FrameLabel->{"Nbr of times a SNV was reported","fraction of SNVs","206 AML SNV locations"},ImageSize->400,ImageMargins->5]

Get patient IDs from library IDs

snvs[[1]]

libIdsAndPatientIdsForSNVs=Table[
{snvs[[i,9]],
Which[
MemberQ[rnaseqLibIds,snvs[[i,9]]]==True,snvs[[i,9]]/.rnaSeqLibId2PatientIDRules,
MemberQ[normalWashULibIDs,snvs[[i,9]]]==True,snvs[[i,9]]/.normalWashULib2PatientRules,
MemberQ[tumorWashULibIDs,snvs[[i,9]]]==True,snvs[[i,9]]/.tumorWashULib2PatientRules,
True,"no rule found"
]}
,{i,Length[snvs]}];
Length[libIdsAndPatientIdsForSNVs]
libIdsAndPatientIdsForSNVs[[1;;3]]//TableForm

uniqPatientIds=Union[libIdsAndPatientIdsForSNVs[[All,2]]]

For a single SNV

s=sortedTallyIds[[1]]
s1=s[[1]]
selS1=Select[snvs,(#[[2]]==s1)&];
librariesS1=Tally[selS1[[All,9]]]

ids=librariesS1[[All,1]]

StringMatchQ["A00002",RegularExpression["\\bA\\d+"]]
StringMatchQ["A00002_N",RegularExpression["\\bA\\d+_N"]]
StringMatchQ["A00002_N",RegularExpression["\\bA\\d+_T"]]
StringMatchQ["A00002_T",RegularExpression["\\bA\\d+_T"]]
(*rnaSeqRegEx="\\bA\\d+"
genomicNormalRegEx="\\bA\\d+_N"
genomicTumorRegEx="\\bA\\d+_T"
libId="A00002"
rnaSeqLibraryQ[libId_]=(StringMatchQ[libId,RegularExpression["\\bA\\d+_N"]])*)

Table[librariesS1[[i,1]],{i,5}]
(*Table[StringMatchQ[librariesS1[[i,1]],RegularExpression["\\bA\\d+"]],{i,Length[librariesS1]}]*)
{Count[Table[StringMatchQ[librariesS1[[i,1]],RegularExpression["\\bA\\d+"]],{i,Length[librariesS1]}],True],"RNA-seq library IDs"}
{Count[Table[StringMatchQ[librariesS1[[i,1]],RegularExpression["\\bA\\d+_N"]],{i,Length[librariesS1]}],True],"normal genomic library IDs"}
{Count[Table[StringMatchQ[librariesS1[[i,1]],RegularExpression["\\bA\\d+_T"]],{i,Length[librariesS1]}],True],"tumor genomic library IDs"}