In this sheet we will load RNAsequencing data for multipl DLBCL cell lines and see what loading that data does to the model.
First let's load the RNA seq data. This has been downloaded from https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103934. Raw counts were normalized to log10-Counts Per Million (logCPM) by dividing each column by the total sum of its counts, multiplying it by 10^6, followed by the application of a log10-transform.
Make sure you have a julia kernel selected that has many threads. Here is how many threads you have available:
import Base.Threads
Threads.nthreads()
12
#packages we need
using DifferentialEquations
using Plots
using CSV
using Distributions
using Random
using DataFrames
using JLD2
using FileIO
using StatsPlots
using Plots.PlotMeasures
using Statistics
plotly()
colorArray=palette(:seaborn_colorblind)
┌ Warning: PlotlyBase 0.5.4 is not compatible with this version of Plots. The declared compatibility is 0.7. └ @ Plots /home/simon/.julia/packages/Plots/yJrrq/src/Plots.jl:22
geneExpressionFileName="RNAseqData/cellLineRNAseqlog10CPM_justDLBCL.csv"
cellLineDF = DataFrame(CSV.File(geneExpressionFileName))
#genesOfInterest=["REL","RELA","RELB","NFKB1","NFKB2","MAP3K14"]
genesOfInterest=["REL","RELA","RELB","NFKB1","NFKB2"]
cellLineDF=filter(row -> (row."Sample Name" in genesOfInterest), cellLineDF)
display(cellLineDF)
Sample Name | DB | DOHH2 | FARAGE | HBL1 | KARPAS422 | OCILY1 | OCILY18 | OCILY19 | |
---|---|---|---|---|---|---|---|---|---|
String15 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | NFKB1 | 2.64493 | 2.91356 | 3.14002 | 2.84536 | 1.98816 | 2.39363 | 2.61325 | 2.55143 |
2 | NFKB2 | 2.15054 | 2.86502 | 3.12984 | 2.63469 | 1.5011 | 2.03163 | 2.15332 | 2.17576 |
3 | REL | 3.0735 | 2.60196 | 3.40004 | 2.85277 | 2.53454 | 2.75115 | 3.10439 | 2.67341 |
4 | RELA | 1.6919 | 1.88115 | 1.71495 | 1.79965 | 1.63255 | 1.48239 | 1.87405 | 1.77258 |
5 | RELB | 2.01452 | 2.14952 | 2.6369 | 2.20556 | 1.00657 | 1.93104 | 1.803 | 1.58667 |
Let's standardize the expression levels for each line to see whether they're over or under expressed.
for rowIndex in 1:size(cellLineDF,1)
meanOfRow=mean(cellLineDF[rowIndex,2:end])
stdOfROw=std(cellLineDF[rowIndex,2:end])
for colIndex in 2:size(cellLineDF,2)
cellLineDF[rowIndex,colIndex]= (cellLineDF[rowIndex,colIndex].-meanOfRow)./stdOfROw
end
end
display(cellLineDF)
Sample Name | DB | DOHH2 | FARAGE | HBL1 | KARPAS422 | OCILY1 | |
---|---|---|---|---|---|---|---|
String15 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | NFKB1 | 0.0579935 | 1.16152 | 2.09178 | 0.881351 | -2.64002 | -0.97436 |
2 | NFKB2 | -0.672191 | 1.10121 | 1.75851 | 0.5295 | -2.28414 | -0.967334 |
3 | REL | 0.724368 | -0.833706 | 1.80335 | -0.00497959 | -1.05648 | -0.340773 |
4 | RELA | -0.217324 | 0.990217 | -0.0702326 | 0.470189 | -0.596045 | -1.55418 |
5 | RELB | -0.298829 | 0.0666828 | 1.38629 | 0.218432 | -3.02786 | -0.524836 |
Let's plot some overall histograms and highlight a few of our favourite lines.
linesToHighlight=["U2932","RIVA","RCK8","DB","OCILY1","OCILY19","OCILY7","OCILY8"]
pairsToPlot=[[3,4],[4,5]]
histList=repeat([plot(1)], length(pairsToPlot))
# for i in 1:length(genesOfInterest)
# thisHist=histogram(Vector(cellLineDF[i,2:end]),title=cellLineDF[i,1],legend=false,bins=:rice,color=colorArray[1],linewidth=2)
# for lineIndex in 1:length(linesToHighlight)
# thisLineIndex=findall(names(cellLineDF).==linesToHighlight[lineIndex])[1]
# plot!(thisHist,[cellLineDF[i,thisLineIndex],cellLineDF[i,thisLineIndex]],[0,20],linewidth=5,color=colorArray[lineIndex],hover=linesToHighlight[lineIndex],dpi=300,size=(500,500))
# end
# plotList[i]=thisHist
# end
for pairIndex in 1:length(pairsToPlot)
thisHist=scatter(Vector(cellLineDF[pairsToPlot[pairIndex][1],2:end]),Vector(cellLineDF[pairsToPlot[pairIndex][2],2:end]),title=cellLineDF[pairsToPlot[pairIndex][1],1]*" vs "*cellLineDF[pairsToPlot[pairIndex][2],1],xlabel=cellLineDF[pairsToPlot[pairIndex][1],1],ylabel=cellLineDF[pairsToPlot[pairIndex][2],1],label=false,color=:black,hover = names(cellLineDF)[2:end],markerstrokewidth=0,markersize=10,dpi=300,size=(500,500))
histList[pairIndex]=thisHist
for lineIndex in 1:length(linesToHighlight)
thisLineIndex=findall(names(cellLineDF).==linesToHighlight[lineIndex])[1]
if pairIndex==1
scatter!(Vector([cellLineDF[pairsToPlot[pairIndex][1],thisLineIndex]]),Vector([cellLineDF[pairsToPlot[pairIndex][2],thisLineIndex]]),color=colorArray[lineIndex],markersize=15,label=linesToHighlight[lineIndex],markerstrokewidth=0)
else
scatter!(Vector([cellLineDF[pairsToPlot[pairIndex][1],thisLineIndex]]),Vector([cellLineDF[pairsToPlot[pairIndex][2],thisLineIndex]]),color=colorArray[lineIndex],markersize=15,label=false,markerstrokewidth=0)
end
end
end
#display(plot(plotList...,size=(1500,800)))
plot(histList...,left_margin = 5mm,bottom_margin=10mm,size=(750,500))
cellLineDF
Sample Name | DB | DOHH2 | FARAGE | HBL1 | KARPAS422 | OCILY1 | |
---|---|---|---|---|---|---|---|
String15 | Float64 | Float64 | Float64 | Float64 | Float64 | Float64 | |
1 | NFKB1 | 0.0579935 | 1.16152 | 2.09178 | 0.881351 | -2.64002 | -0.97436 |
2 | NFKB2 | -0.672191 | 1.10121 | 1.75851 | 0.5295 | -2.28414 | -0.967334 |
3 | REL | 0.724368 | -0.833706 | 1.80335 | -0.00497959 | -1.05648 | -0.340773 |
4 | RELA | -0.217324 | 0.990217 | -0.0702326 | 0.470189 | -0.596045 | -1.55418 |
5 | RELB | -0.298829 | 0.0666828 | 1.38629 | 0.218432 | -3.02786 | -0.524836 |
#linesToHighlight=["U2932","SUDHL10","TMD8","RIVA","RCK8"]
cooFileName="RNAseqData/cellLineRNAseq_justDLBCL_COO.csv"
cooDF = DataFrame(CSV.File(cooFileName))
colorOfPoints=Array(cooDF[1,2:end])
colorArray=zeros(size(colorOfPoints))
colorArray[findall(x->x=="GCB",colorOfPoints)].=1
colorArray[findall(x->x=="ABC",colorOfPoints)].=2
CSV.write("cellLineDF.csv",cellLineDF)
pairsToPlot=[[1,2]]
plotList = repeat([plot(1)], length(genesOfInterest))
histList=repeat([plot(1)], length(pairsToPlot))
for pairIndex in 1:length(pairsToPlot)
thisHist=scatter(Vector(cellLineDF[pairsToPlot[pairIndex][1],2:end]),Vector(cellLineDF[pairsToPlot[1][2],2:end]),title=cellLineDF[pairsToPlot[pairIndex][1],1]*" vs "*cellLineDF[pairsToPlot[pairIndex][2],1],xlabel=cellLineDF[pairsToPlot[pairIndex][1],1],ylabel=cellLineDF[pairsToPlot[pairIndex][2],1],label=false,group=colorArray[1:end],hover = names(cooDF)[2:end],markerstrokewidth=0,markersize=15,dpi=300,size=(500,500))
histList[pairIndex]=thisHist
end
plot(histList...,left_margin = 5mm,bottom_margin=10mm,size=(600,600))
cellLineDF[:,1]
5-element Vector{String15}: "NFKB1" "NFKB2" "REL" "RELA" "RELB"
cellLineDF.AVG = repeat([0.0],length(genesOfInterest))
linesToHighlight=linesToHighlight[1:2]
#"REL","RELA","RELB","NFKB1","NFKB2","MAP3K14"
#variableList=["basal_cRelSynth","k1_RelASynth","basal_RelBSynth","basal_p50Synth","kmax_p100Synth","k1_NIKSynth"]
variableList=["basal_p50Synth","kmax_p100Synth","basal_cRelSynth","k1_RelASynth","basal_RelBSynth"]
conditions=linesToHighlight
numberOfLines=length(conditions)
paramsToChange=repeat([variableList],numberOfLines)
modifyAmount=repeat([[1.0]],numberOfLines)
for i in 1:numberOfLines
thisLineIndex=findall(names(cellLineDF).==linesToHighlight[i])[1]
scalingVals=Vector(cellLineDF[:,thisLineIndex])
scalingVals=10 .^scalingVals
modifyAmount[i]=scalingVals
end
first_cell=1
last_cell=25
#set up where CSV2Julia is
#from: https://github.com/SiFTW/CSV2JuliaDiffEq
locationOfCSV2Julia="CSV2Julia/csv2model-multiscale.py"
#identify the three CSV sheets that describe the model
reactionsFile="moduleDefinitionFiles/updated_NFkB_reactions.csv"
parametersFile="moduleDefinitionFiles/parameters.csv"
rateLawsFile="moduleDefinitionFiles/rateLaws.csv"
generatedCSVLocation="generatedCSVs/"
distributedModelFilesLocation="distributedModelFiles/"
mkpath(generatedCSVLocation)
mkpath(distributedModelFilesLocation)
colorArray=palette(:seaborn_colorblind)
totalIKK=140
total_WT_IKK=totalIKK*1.0
basalIKK=totalIKK/100
IKKMultiplier=1
maxTimeSS=100000.0
maxTimeTC=60*8
preCV=0.11;
#preCV=0.0
function ikkDefault(t,maxTime)
#explicitely defined input
IKKparamVals=[0,85,100,85,45,30,25,23,18,15,12,10,10]./100
paramTime=[0,5,10,15,18,20,300,400,450,480,500,520,max(maxTime,2880)]
#get the value after the current time point
indexGTTime=min(searchsortedfirst(paramTime,t),length(paramTime))
#get the one before
indexLTTime=max(1,indexGTTime-1)
timeGTt=paramTime[indexGTTime]
timeLTt=paramTime[indexLTTime]
valGTt=IKKparamVals[indexGTTime]
valLTt=IKKparamVals[indexLTTime]
timeDiff=timeGTt-timeLTt
#if the time before and time after are different do a basic linear interpolation between the two.
if timeDiff>0
valDiff=valGTt-valLTt
gradient=valDiff/timeDiff
timeStep=t-timeLTt
return valLTt+(gradient*timeStep)
else
return valLTt
end
end
function NIKDefault(t,maxTime)
#explicitely defined input
NIKparamVals=[1,50,100,180,200,100,50]./100
paramTime=[0,60,120,300,400,800,max(maxTime,2880)]
#get the value after the current time point
indexGTTime=min(searchsortedfirst(paramTime,t),length(paramTime))
#get the one before
indexLTTime=max(1,indexGTTime-1)
timeGTt=paramTime[indexGTTime]
timeLTt=paramTime[indexLTTime]
valGTt=NIKparamVals[indexGTTime]
valLTt=NIKparamVals[indexLTTime]
timeDiff=timeGTt-timeLTt
#if the time before and time after are different do a basic linear interpolation between the two.
if timeDiff>0
valDiff=valGTt-valLTt
gradient=valDiff/timeDiff
timeStep=t-timeLTt
return valLTt+(gradient*timeStep)
else
return valLTt
end
end
#the IKK function is just maps to basal IKK during SS and basal+the IKK curve during the time course.
ikkSS=t->basalIKK
ikkTC=t->basalIKK+(ikkDefault(t, maxTimeTC)*IKKMultiplier*totalIKK);
ikkTCWT=t->basalIKK+(ikkDefault(t, maxTimeTC)*IKKMultiplier*total_WT_IKK);
ikkCurveForLine=[ikkTCWT]
ikkFunc=ikkTCWT
NIKFuncSS=t->1
NIKFuncTC=t->0
ikkSSHigh=t->basalIKKHigh
p=plot([1:1:maxTimeTC],ikkTC.(1:1:maxTimeTC),label="smooth fit",title="input curve");
xlabel!("time (h)")
ylabel!("active IKK")
The below code first creates the ODE file.
parametersDF = DataFrame(CSV.File(parametersFile,types=Dict(:parameter=>String, :value=>String, :distribute=>Int64)))
originalParams=deepcopy(parametersDF)
thisModelName="odeModel.jl"
thisParamFile=parametersFile
arguments=[reactionsFile, thisParamFile, rateLawsFile,thisModelName]
cmd=`python3 $locationOfCSV2Julia $arguments param`
#lets run csv2julia (requires python to be installed)
run(cmd)
include(thisModelName)
include("variableNames.jl")
include("scanIncludes.jl")
mv(thisModelName,"distributedModelFiles/"*thisModelName, force=true)
println("Model generated for all conditions")
param Running CSV2JuliaDiffEq with parameters dynamically determined by a variable, re-run with the 5th argument set to 'scan' or 'inline' Opening moduleDefinitionFiles/rateLaws.csv as rate law file Opening moduleDefinitionFiles/parameters.csv as parameters file Opening moduleDefinitionFiles/updated_NFkB_reactions.csv as reactions file parameters can now be searched in parameterNameList by name. example to modify k_binding 1.5 fold higher: indexOfParam=findfirst(x->"k_binding"==x,parameterNameList) paramVals[indexOfParam]=paramVals[indexOfParam]*1.5 Model generated for all conditions
The below function will run the number of cells you specify, with the conditions you specify and the input curves specified. It will use multiple threads to do so and save the results in a different folder for each condition.
function runSimulationNew(first_cell, last_cell, conditions,folder,IKKSSArray,IKKTCArray,NIKSSArray,NIKTCArray)
mkpath(folder)
#now lets loop through and solve the cell
TCLength=1000*60
maximumAttemptsAtSS=10
originalParams=copy(paramVals)
ikkIndex=findfirst(x -> x=="ikk1_ikkactivity", parameterNameList)
NIKIndex=findfirst(x -> x=="nik_deg_mod", parameterNameList)
thisDist=TruncatedNormal(1.0, preCV,0,Inf)
Random.seed!(123)
allParams=[]
allParams=Array{Any}(undef, size(parametersDF,1), last_cell)
for cellIndex in first_cell:last_cell
thisCellsParamVals=copy(originalParams)
for j in 1:size(parametersDF,1)
if parametersDF[j,3]==1
x = rand(thisDist, 1)
thisCellsParamVals[j]=thisCellsParamVals[j].*x[1]
end
end
#println(thisCellsParamVals)
allParams[:,cellIndex]=thisCellsParamVals
end
df = DataFrame(allParams,:auto)
#add the variable names and save to a file
insertcols!(df, 1, :names=>parameterNameList)
CSV.write(folder*"/allParams_runSimulationNew.csv",df);
#println(size(allParams))
#println(allParams[1])
allParamsOriginal=copy(allParams)
for condIndex in 1:length(conditions)
allParams=copy(allParamsOriginal)
thisCondition = conditions[condIndex]
#TODO: consider making a condition scaling factor array here and just multiplying all prameters by it every time.
println("Starting condition: "*thisCondition)
odeName="odeModel"
myFun=getfield(Main,Symbol(odeName))
#define the function and the initial conditions
f=ODEFunction(myFun,syms=Symbol.(syms))
y0=zeros(length(syms))
paramsListInThisCondition=paramsToChange[condIndex]
modifyListInThisCondition=modifyAmount[condIndex]
for thisParamIndex in 1:length(paramsListInThisCondition)
thisParam=paramsListInThisCondition[thisParamIndex]
thisParamsIndexInParamList=findfirst(x->x==thisParam,parameterNameList)
allParams[thisParamsIndexInParamList,:]=allParams[thisParamsIndexInParamList,:].*modifyListInThisCondition[thisParamIndex]
end
df = DataFrame(allParams,:auto)
#add the variable names and save to a file
insertcols!(df, 1, :names=>parameterNameList)
CSV.write(folder*"/allParams_runSimulationNew_"*thisCondition*".csv",df);
Threads.@threads for i in first_cell:last_cell
#figure out the name of this cell's ode file
#DISTRIBUTE PARAMS
println("starting cell: "*string(i))
thisCellsParamVals=copy(allParams[:,i])
thisCellsParamVals[ikkIndex]=IKKSSArray[condIndex]
thisCellsParamVals[NIKIndex]=NIKSSArray[condIndex]
# #now write this condition's CSV file to a folder of cells
# CSV.write(generatedCSVLocation*"/parameters_"*string(conditions[condIndex])*".csv", thisCondParamFile)
prob=ODEProblem(f,y0,(0.0,maxTimeSS),thisCellsParamVals)
solss=solve(prob,saveat=100.0,progress = true)
println("Steady state found for cell: "*string(i))
#dynamic phase, use SS solution as initial conditions
y0=vec(convert(Array, solss[:,end]))
y0[y0.<0].=0
thisCellsParamVals[ikkIndex]=IKKTCArray[condIndex]
thisCellsParamVals[NIKIndex]=NIKTCArray[condIndex]
#CSV.write(generatedCSVLocation*"/parameters_"*string(conditions[condIndex])*"_cell_"*string(i)*".csv", DataFrame(thisCellsParamVals,:auto))
try
f=ODEFunction(myFun,syms=syms)
prob=ODEProblem(f,y0,(0.0,maxTimeTC),thisCellsParamVals)
println("Solving equations for dynamic time course for cell:"*string(i))
sol=solve(prob, abstol=1e-5,reltol=1e-3, saveat=1.0)
#save("outputs/sol_"*thisCondition*"_cell_"*string(i)*".jld2", "solution", sol)
df = DataFrame(Float64.(sol),:auto)
#add the variable names and save to a file
insertcols!(df, 1, :names=>syms)
#CSV.write("outputs/sol_"*thisCondition*"_cell_"*string(i)*".csv",Tables.columntable(df));
CSV.write(folder*"/sol_"*thisCondition*"_cell_"*string(i)*".csv",df);
catch e
println("error:")
println(e)
end
end
println("all cells done in condition: "*thisCondition)
end
end
runSimulationNew (generic function with 1 method)
The below code is used to plot some graphs. It will plot the species you specify and compare between each condition. It will plot both steady state values (as bar graphs on the left), and time courses (as line graphs on the right). Then
function plotAllSpecies(speciesToPlot,conditionsToPlot,colorArray,first_cell,last_cell,folder,hoursToPlot)
for species in speciesToPlot
#thisPlot=plot(title=species)
thisPlotStd=plot(title=species*" avg")
boxPlotAll=plot(title=species*" ss")
BoxPlotAvg=plot(title=species*"avg ss")
conditionIndex=1
meansOfAllConditions=zeros(1,length(conditionsToPlot))
stdOfAllConditions=zeros(1,length(conditionsToPlot))
lengthOfTC=0
for condition in conditionsToPlot
lengthOfTC=size(DataFrame(CSV.File(folder*"/sol_"*condition*"_cell_1.csv")),2)-1
conditionArray=zeros(last_cell,lengthOfTC)
lineColor=colorArray[conditionIndex]
virtExpFlag=false
for i in first_cell:last_cell
thisCellData=DataFrame(CSV.File(folder*"/sol_"*condition*"_cell_"*string(i)*".csv"))
if !("names" in names(thisCellData))
insertcols!(thisCellData, 1, :names=>syms)
end
allNoneFloats=findall(eltype.(eachcol(thisCellData)).!=Float64)
if length(allNoneFloats)>1
for index in allNoneFloats[2:end]
thisCellData[!,index]=parse.(Float64,thisCellData[:,index])
end
end
thisTC=zeros(1,size(thisCellData,2)-1)
if endswith(species,"*")
virtExpFlag=true
speciesShort=species[1:end-1]
speciesIDs=intersect(findall( x ->occursin(speciesShort,x),syms),findall(x->!startswith(x,"t"),syms))
speciesNames=String.(syms[speciesIDs])
# println("For species: "*species*" printing: ")
# println(speciesNames)
for thisName in speciesNames
thisSpeciesTC=convert(Matrix, thisCellData[thisCellData[!,:names].==thisName,:])[2:end]
thisTC=thisTC.+thisSpeciesTC'
end
else
thisTC=convert(Matrix, thisCellData[thisCellData[!,:names].==species,:])[2:end]
end
conditionArray[i,:]=thisTC[1:lengthOfTC]
end
df = DataFrame(Float64.(conditionArray),:auto)
#add the variable names and save to a file
#CSV.write("outputs/sol_"*thisCondition*"_cell_"*string(i)*".csv",Tables.columntable(df));
CSV.write("outputs/allTCs_"*species*"_cell.csv",df);
meanOfCondition=mean(conditionArray, dims=1)
stdOfCondition=std(conditionArray, dims=1)
# println(meanOfCondition)
# println(stdOfCondition)
plot!(thisPlotStd,meanOfCondition',grid=false,color=lineColor,ribbon=stdOfCondition',fillalpha=.5,label=condition,linewidth=5)
meansOfAllConditions[conditionIndex]=meanOfCondition[1]
stdOfAllConditions[conditionIndex]=stdOfCondition[1]
conditionIndex+=1
end
conditionIndex=1
#plot!(boxPlotAll,conditionsToPlot, meansOfAllConditions;, c=colorArray, yerr = stdOfAllConditions', label = "",xrotation = 90,seriestype = :scatter,fillcolor=:match)
for condition in conditionsToPlot
plot!(boxPlotAll,[conditionIndex], [meansOfAllConditions[conditionIndex]], c=colorArray[conditionIndex], yerr = stdOfAllConditions[conditionIndex], label = false,xrotation = 90,seriestype = :scatter,fillcolor=:match,markersize=20,markerstrokewidth=5)
conditionIndex+=1
end
plot!(boxPlotAll,xticks = (1:length(conditions), conditions),xlim=(0,length(conditionsToPlot)+1),ylim=(0,maximum([maximum(meansOfAllConditions)+maximum(stdOfAllConditions),10])),dpi=300,size=(800,1000),xtickfontsize=18,ytickfontsize=18)
plot!(thisPlotStd,xticks=(collect(0:30:lengthOfTC),collect(0:0.5:lengthOfTC/60)),ylim=(0,maximum(meansOfAllConditions)+maximum(stdOfAllConditions)),xlim=(0,hoursToPlot*60),dpi=300,size=(1500,1000),xtickfontsize=18,ytickfontsize=18)
#display(plot(boxPlotAll,thisPlot))
display(plot(boxPlotAll,thisPlotStd,layout = grid(1, 2, widths=[0.4 ,0.6])))
end
end
plotAllSpecies (generic function with 1 method)
function plotAllSpeciesNormalised(speciesToPlot,conditionsToPlot,colorArray,first_cell,last_cell,folder,hoursToPlot,ymax)
for species in speciesToPlot
#thisPlot=plot(title=species)
# thisPlotStd=plot(title=species*" avg")
# boxPlotAll=plot(title=species*" ss")
# BoxPlotAvg=plot(title=species*"avg ss")
thisPlotStd=plot()
boxPlotAll=plot()
BoxPlotAvg=plot()
conditionIndex=1
meansOfAllConditions=zeros(1,length(conditionsToPlot))
stdOfAllConditions=zeros(1,length(conditionsToPlot))
lengthOfTC=0
for condition in conditionsToPlot
lengthOfTC=size(DataFrame(CSV.File(folder*"/sol_"*condition*"_cell_1.csv")),2)-1
conditionArray=zeros(last_cell,lengthOfTC)
lineColor=colorArray[conditionIndex]
virtExpFlag=false
for i in first_cell:last_cell
thisCellData=DataFrame(CSV.File(folder*"/sol_"*condition*"_cell_"*string(i)*".csv"))
if !("names" in names(thisCellData))
insertcols!(thisCellData, 1, :names=>syms)
end
allNoneFloats=findall(eltype.(eachcol(thisCellData)).!=Float64)
if length(allNoneFloats)>1
for index in allNoneFloats[2:end]
thisCellData[!,index]=parse.(Float64,thisCellData[:,index])
end
end
thisTC=zeros(1,size(thisCellData,2)-1)
if endswith(species,"*")
virtExpFlag=true
speciesShort=species[1:end-1]
speciesIDs=intersect(findall( x ->occursin(speciesShort,x),syms),findall(x->!startswith(x,"t"),syms))
speciesNames=String.(syms[speciesIDs])
# println("For species: "*species*" printing: ")
# println(speciesNames)
for thisName in speciesNames
thisSpeciesTC=convert(Matrix, thisCellData[thisCellData[!,:names].==thisName,:])[2:end]
thisTC=thisTC.+thisSpeciesTC'
end
else
thisTC=convert(Matrix, thisCellData[thisCellData[!,:names].==species,:])[2:end]
end
conditionArray[i,:]=thisTC[1:lengthOfTC]
end
df = DataFrame(Float64.(conditionArray),:auto)
#add the variable names and save to a file
#CSV.write("outputs/sol_"*thisCondition*"_cell_"*string(i)*".csv",Tables.columntable(df));
CSV.write("outputs/allTCs_"*species*"_cell.csv",df);
meanOfCondition=mean(conditionArray, dims=1)
stdOfCondition=std(conditionArray, dims=1)
# println(meanOfCondition)
# println(stdOfCondition)
plot!(thisPlotStd,meanOfCondition'./meanOfCondition[1],grid=false,color=lineColor,ribbon=stdOfCondition'./meanOfCondition[1],fillalpha=.5,label=condition,linewidth=5)
meansOfAllConditions[conditionIndex]=meanOfCondition[1]
stdOfAllConditions[conditionIndex]=stdOfCondition[1]
conditionIndex+=1
end
conditionIndex=1
#plot!(boxPlotAll,conditionsToPlot, meansOfAllConditions;, c=colorArray, yerr = stdOfAllConditions', label = "",xrotation = 90,seriestype = :scatter,fillcolor=:match)
for condition in conditionsToPlot
plot!(boxPlotAll,[conditionIndex], [meansOfAllConditions[conditionIndex]], c=colorArray[conditionIndex], yerr = stdOfAllConditions[conditionIndex], label = false,xrotation = 90,seriestype = :scatter,fillcolor=:match,markersize=8,markerstrokewidth=3)
conditionIndex+=1
end
plot!(boxPlotAll,xticks = (1:length(conditions), conditions),xlim=(0,length(conditionsToPlot)+1),ylim=(0,maximum([maximum(meansOfAllConditions)+maximum(stdOfAllConditions).+1,10])),dpi=300,size=(200,500),xtickfontsize=18,ytickfontsize=18)
plot!(thisPlotStd,xticks=(collect(0:30:lengthOfTC),collect(0:0.5:lengthOfTC/60)),ylim=(0,ymax),xlim=(0,hoursToPlot*60),dpi=300,size=(750,500),xtickfontsize=18,ytickfontsize=18)
#display(plot(boxPlotAll,thisPlot))
display(plot(boxPlotAll,thisPlotStd,layout = grid(1, 2, widths=[0.2 ,0.8])))
end
end
plotAllSpeciesNormalised (generic function with 1 method)
function plotAllSpeciesInOneLineNormalised(speciesToPlot,linesToPlot,colorArray,first_cell,last_cell,folder,hoursToPlot,ymax)
for lines in linesToPlot
#thisPlot=plot(title=species)
thisPlotStd=plot(title=lines*" avg")
boxPlotAll=plot(title=lines*" ss")
specieIndex=1
meansOfAllSpecies=zeros(1,length(speciesToPlot))
stdOfAllSpecies=zeros(1,length(speciesToPlot))
lengthOfTC=0
for species in speciesToPlot
lengthOfTC=size(DataFrame(CSV.File(folder*"/sol_"*lines*"_cell_1.csv")),2)-1
cellLineArray=zeros(last_cell,lengthOfTC)
lineColor=colorArray[specieIndex]
for i in first_cell:last_cell
thisCellData=DataFrame(CSV.File(folder*"/sol_"*lines*"_cell_"*string(i)*".csv"))
if !("names" in names(thisCellData))
insertcols!(thisCellData, 1, :names=>syms)
end
allNoneFloats=findall(eltype.(eachcol(thisCellData)).!=Float64)
if length(allNoneFloats)>1
for index in allNoneFloats[2:end]
thisCellData[!,index]=parse.(Float64,thisCellData[:,index])
end
end
thisTC=zeros(1,size(thisCellData,2)-1)
if endswith(species,"*")
virtExpFlag=true
speciesShort=species[1:end-1]
speciesIDs=intersect(findall( x ->occursin(speciesShort,x),syms),findall(x->!startswith(x,"t"),syms))
speciesNames=String.(syms[speciesIDs])
# println("For species: "*species*" printing: ")
# println(speciesNames)
for thisName in speciesNames
thisSpeciesTC=convert(Matrix, thisCellData[thisCellData[!,:names].==thisName,:])[2:end]
thisTC=thisTC.+thisSpeciesTC'
end
else
thisTC=convert(Matrix, thisCellData[thisCellData[!,:names].==species,:])[2:end]
end
cellLineArray[i,:]=thisTC[1:lengthOfTC]
end
df = DataFrame(Float64.(cellLineArray),:auto)
#add the variable names and save to a file
#CSV.write("outputs/sol_"*thisCondition*"_cell_"*string(i)*".csv",Tables.columntable(df));
CSV.write("outputs/allTCs_"*species*"_"*lines*"_cell.csv",df);
meanOfSpecies=mean(cellLineArray, dims=1)
stdOfSpecies=std(cellLineArray, dims=1)
plot!(thisPlotStd,meanOfSpecies'./meanOfSpecies[1],grid=false,color=lineColor,ribbon=stdOfSpecies'./meanOfSpecies[1],fillalpha=.5,label=species,linewidth=5)
meansOfAllSpecies[specieIndex]=meanOfSpecies[1]
stdOfAllSpecies[specieIndex]=stdOfSpecies[1]
specieIndex+=1
end
#plot!(boxPlotAll,conditionsToPlot, meansOfAllConditions;, c=colorArray, yerr = stdOfAllConditions', label = "",xrotation = 90,seriestype = :scatter,fillcolor=:match)
specieIndex=1
for species in speciesToPlot
plot!(boxPlotAll,[specieIndex], [meansOfAllSpecies[specieIndex]], c=colorArray[specieIndex], yerr = stdOfAllSpecies[specieIndex], label = false,xrotation = 90,seriestype = :scatter,fillcolor=:match,markersize=20,markerstrokewidth=5)
specieIndex+=1
end
plot!(thisPlotStd,xticks=(collect(0:30:lengthOfTC),collect(0:0.5:lengthOfTC/60)),xlim=(0,hoursToPlot*60),ylim=(0,ymax),dpi=300,size=(1500,1000),xtickfontsize=18,ytickfontsize=18)
plot!(boxPlotAll,xticks = (1:length(speciesToPlot), speciesToPlot),xlim=(0,length(speciesToPlot)+1),ylim=(0,maximum([maximum(meansOfAllSpecies)+maximum(stdOfAllSpecies),10])),dpi=300,size=(500,1000),xtickfontsize=18,ytickfontsize=18)
#display(plot(boxPlotAll,thisPlot))
display(plot(boxPlotAll,thisPlotStd,layout = grid(1, 2, widths=[0.25 ,0.75],dpi=300)))
#savefig(lines*"_selectSpecies.png")
end
end
plotAllSpeciesInOneLineNormalised (generic function with 1 method)
function plotAllSpeciesInOneLine(speciesToPlot,linesToPlot,colorArray,first_cell,last_cell,folder,hoursToPlot)
for lines in linesToPlot
#thisPlot=plot(title=species)
thisPlotStd=plot(title=lines*" avg")
boxPlotAll=plot(title=lines*" ss")
specieIndex=1
meansOfAllSpecies=zeros(1,length(speciesToPlot))
stdOfAllSpecies=zeros(1,length(speciesToPlot))
lengthOfTC=0
for species in speciesToPlot
lengthOfTC=size(DataFrame(CSV.File(folder*"/sol_"*lines*"_cell_1.csv")),2)-1
cellLineArray=zeros(last_cell,lengthOfTC)
lineColor=colorArray[specieIndex]
for i in first_cell:last_cell
thisCellData=DataFrame(CSV.File(folder*"/sol_"*lines*"_cell_"*string(i)*".csv"))
if !("names" in names(thisCellData))
insertcols!(thisCellData, 1, :names=>syms)
end
allNoneFloats=findall(eltype.(eachcol(thisCellData)).!=Float64)
if length(allNoneFloats)>1
for index in allNoneFloats[2:end]
thisCellData[!,index]=parse.(Float64,thisCellData[:,index])
end
end
thisTC=zeros(1,size(thisCellData,2)-1)
if endswith(species,"*")
virtExpFlag=true
speciesShort=species[1:end-1]
speciesIDs=intersect(findall( x ->occursin(speciesShort,x),syms),findall(x->!startswith(x,"t"),syms))
speciesNames=String.(syms[speciesIDs])
# println("For species: "*species*" printing: ")
# println(speciesNames)
for thisName in speciesNames
thisSpeciesTC=convert(Matrix, thisCellData[thisCellData[!,:names].==thisName,:])[2:end]
thisTC=thisTC.+thisSpeciesTC'
end
else
thisTC=convert(Matrix, thisCellData[thisCellData[!,:names].==species,:])[2:end]
end
cellLineArray[i,:]=thisTC[1:lengthOfTC]
end
df = DataFrame(Float64.(cellLineArray),:auto)
#add the variable names and save to a file
#CSV.write("outputs/sol_"*thisCondition*"_cell_"*string(i)*".csv",Tables.columntable(df));
CSV.write("outputs/allTCs_"*species*"_"*lines*"_cell.csv",df);
meanOfSpecies=mean(cellLineArray, dims=1)
stdOfSpecies=std(cellLineArray, dims=1)
plot!(thisPlotStd,meanOfSpecies',grid=false,color=lineColor,ribbon=stdOfSpecies',fillalpha=.5,label=species,linewidth=5)
meansOfAllSpecies[specieIndex]=meanOfSpecies[1]
stdOfAllSpecies[specieIndex]=stdOfSpecies[1]
specieIndex+=1
end
#plot!(boxPlotAll,conditionsToPlot, meansOfAllConditions;, c=colorArray, yerr = stdOfAllConditions', label = "",xrotation = 90,seriestype = :scatter,fillcolor=:match)
specieIndex=1
for species in speciesToPlot
plot!(boxPlotAll,[specieIndex], [meansOfAllSpecies[specieIndex]], c=colorArray[specieIndex], yerr = stdOfAllSpecies[specieIndex], label = false,xrotation = 90,seriestype = :scatter,fillcolor=:match,markersize=20,markerstrokewidth=5)
specieIndex+=1
end
plot!(thisPlotStd,xticks=(collect(0:60:lengthOfTC),collect(0:1:lengthOfTC/60)),xlim=(0,hoursToPlot*60),ylim=(0,150),dpi=300,size=(1500,1000),xtickfontsize=18,ytickfontsize=18)
plot!(boxPlotAll,xticks = (1:length(speciesToPlot), speciesToPlot),xlim=(0,length(speciesToPlot)+1),ylim=(0,maximum([maximum(meansOfAllSpecies)+maximum(stdOfAllSpecies),10])),dpi=300,size=(500,1000),xtickfontsize=18,ytickfontsize=18)
#display(plot(boxPlotAll,thisPlot))
display(plot(boxPlotAll,thisPlotStd,layout = grid(1, 2, widths=[0.25 ,0.75],dpi=300)))
#savefig(lines*"_selectSpecies.png")
end
end
plotAllSpeciesInOneLine (generic function with 1 method)
IKKSSArray=repeat([ikkSS],numberOfLines)
IKKTCArray=repeat([ikkTC],numberOfLines)
NIKSSArray=repeat([NIKFuncSS],numberOfLines)
NIKTCArray=repeat([NIKFuncSS],numberOfLines)
println("Summary of conditions being run:")
folder="outputs/cellLineScans"
#conditions=conditions[1:2]
#paramsToChange=paramsToChange[1:2]
#modifyAmount=modifyAmount[1:2]
show(IOContext(stdout, :limit => false), "text/plain", hcat(conditions,paramsToChange,modifyAmount))
Summary of conditions being run: 2×3 Matrix{Any}: "U2932" ["basal_p50Synth", "kmax_p100Synth", "basal_cRelSynth", "k1_RelASynth", "basal_RelBSynth"] [4.51869, 0.837478, 0.747799, 1.10961, 1.80649] "RIVA" ["basal_p50Synth", "kmax_p100Synth", "basal_cRelSynth", "k1_RelASynth", "basal_RelBSynth"] [0.670692, 0.972792, 13.6844, 0.0829159, 2.88281]
If the conditions above look good then run the simulation:
runSimulationNew(first_cell,last_cell,conditions,folder,IKKSSArray,IKKTCArray,NIKSSArray,NIKTCArray)
Starting condition: U2932 starting cell: 1 starting cell: 12 starting cell: 4 starting cell: 8 starting cell: 16 starting cell: 10 starting cell: 24 starting cell: 22 starting cell: 20 starting cell: 14 starting cell: 6 starting cell: 18 Steady state found for cell: 22 Steady state found for cell: 8 Steady state found for cell: 18 Steady state found for cell: 14 Steady state found for cell: 12 Steady state found for cell: 24 Steady state found for cell: 20 Steady state found for cell: 1 Steady state found for cell: 6 Steady state found for cell: 4 Solving equations for dynamic time course for cell:6 Solving equations for dynamic time course for cell:4 Solving equations for dynamic time course for cell:22 Solving equations for dynamic time course for cell:12 Solving equations for dynamic time course for cell:20 Steady state found for cell: 16 Solving equations for dynamic time course for cell:16 Steady state found for cell: 10 Solving equations for dynamic time course for cell:10 Solving equations for dynamic time course for cell:24 Solving equations for dynamic time course for cell:8 Solving equations for dynamic time course for cell:1 Solving equations for dynamic time course for cell:14 Solving equations for dynamic time course for cell:18 starting cell: 11 starting cell: 17 starting cell: 7 starting cell: 13 starting cell: 2 starting cell: 9 starting cell: 25 starting cell: 5 starting cell: 23 starting cell: 21 starting cell: 15 starting cell: 19 Steady state found for cell: 15 Steady state found for cell: 25 Steady state found for cell: 11 Steady state found for cell: 17 Solving equations for dynamic time course for cell:15 Solving equations for dynamic time course for cell:25 Solving equations for dynamic time course for cell:17 Solving equations for dynamic time course for cell:11 Steady state found for cell: 13 Solving equations for dynamic time course for cell:13 Steady state found for cell: 7 Solving equations for dynamic time course for cell:7 Steady state found for cell: 21 Solving equations for dynamic time course for cell:21 Steady state found for cell: 2 Solving equations for dynamic time course for cell:2 Steady state found for cell: 23 Solving equations for dynamic time course for cell:23 Steady state found for cell: 19 Solving equations for dynamic time course for cell:19 Steady state found for cell: 5 Solving equations for dynamic time course for cell:5 Steady state found for cell: 9 Solving equations for dynamic time course for cell:9 starting cell: 3 Steady state found for cell: 3 Solving equations for dynamic time course for cell:3 all cells done in condition: U2932 Starting condition: RIVA starting cell: 1 starting cell: 10 starting cell: 14 starting cell: 6 starting cell: 24 starting cell: 12 starting cell: 16 starting cell: 22 starting cell: 18 starting cell: 20 starting cell: 4 starting cell: 8 Steady state found for cell: 20 Steady state found for cell: 22 Solving equations for dynamic time course for cell:20 Solving equations for dynamic time course for cell:22 Steady state found for cell: 1 Solving equations for dynamic time course for cell:1 starting cell: 21 Steady state found for cell: 4 Solving equations for dynamic time course for cell:4 Steady state found for cell: 10 Solving equations for dynamic time course for cell:10 Steady state found for cell: 16 Solving equations for dynamic time course for cell:16 Steady state found for cell: 8 Solving equations for dynamic time course for cell:8 Steady state found for cell: 18 Solving equations for dynamic time course for cell:18 Steady state found for cell: 14 Solving equations for dynamic time course for cell:14 starting cell: 23 Steady state found for cell: 12 Solving equations for dynamic time course for cell:12 Steady state found for cell: 6 Solving equations for dynamic time course for cell:6 Steady state found for cell: 24 Solving equations for dynamic time course for cell:24 Steady state found for cell: 21 Solving equations for dynamic time course for cell:21 Steady state found for cell: 23 Solving equations for dynamic time course for cell:23 starting cell: 15 starting cell: 19 starting cell: 13 starting cell: 11 starting cell: 9 starting cell: 25 starting cell: 7 starting cell: 2 starting cell: 17 starting cell: 5 Steady state found for cell: 13 Solving equations for dynamic time course for cell:13 Steady state found for cell: 17 Solving equations for dynamic time course for cell:17 Steady state found for cell: 15 Solving equations for dynamic time course for cell:15 Steady state found for cell: 2 Solving equations for dynamic time course for cell:2 Steady state found for cell: 11 Solving equations for dynamic time course for cell:11 Steady state found for cell: 7 Solving equations for dynamic time course for cell:7 Steady state found for cell: 25 Solving equations for dynamic time course for cell:25 Steady state found for cell: 19 Solving equations for dynamic time course for cell:19 Steady state found for cell: 5 Solving equations for dynamic time course for cell:5 Steady state found for cell: 9 Solving equations for dynamic time course for cell:9 starting cell: 3 Steady state found for cell: 3 Solving equations for dynamic time course for cell:3 all cells done in condition: RIVA
colorArray=palette(:seaborn_colorblind)[[1,2]]
folder="outputs/cellLineScans"
linesToHighlight=["U2932","RIVA"]
conditions=linesToHighlight
conditions=conditions[1:end]
speciesToPlot=["cRelnp50n","RelAnp50n"]
plotAllSpeciesNormalised(speciesToPlot,conditions,colorArray,first_cell,last_cell,folder,2,50)