greifflab / immunesim Goto Github PK

Hi,

I am trying immuneSIM with some non-human non-mouse data. I installed the software and can reproduce the code that creates synthetic VDJs and then runs the simulation:

https://immunesim.readthedocs.io/en/latest/parameters.html#introducing-new-germline-gene-frequencies

I can modify the Ds and Js with the sequence from the genome I am interested, and the simulation runs in almost the same time as with the example.

I can modify one of the two V sequences in the example, and leave the frequence of the Vs at 0.5+0.5, and the simulation takes a bit longer but succeeds.

But if I change both Vs for my species of interest, or change only one of the Vs but give it a high frequency of 0.99, then the simulation never ends and uses up 100% of CPU. I left it running for minutes and it didn't complete.

Could it be that it's trying to match the Vs to some other empirical data that is in the code and fails to match it, running in an endless loop?

Which parameter should I switch on/off to be able to run on my own Vs?

Thanks in advance

Hi

I am using immuneSIM to simulate a repertoire with an altered usage of V , D and J genes, but the resulting repertoire V,D,J frequencies do not match the ones I gave in input.

I have altered the frequencies by copying the default "list_germline_genes_allele_01" and substituting the values in the column "frequency".
In the example below I only show the J genes

vdj_list <- list_germline_genes_allele_01
vdj_list$hs$ig$h$J[,"frequency"] <- c( 0.01280559,0.02328289 ,0.08614668, 0.29569267, 0.10942957 ,0.47264261)

sim_repertoire <- immuneSIM(
                number_of_seqs = 500  ,
                vdj_list = vdj_list,
                species = "hs",
                receptor = "ig",
                chain = "h",
                verbose=TRUE,
                shm.mode = "none",
                shm.prob = 7/350,
                freq_update_time=c(250 ),
                vdj_noise = 0,                
                max_cdr3_length = 45)

But when i check the frequencies of J gene, they are quite different from the ones that i set.

> prop.table(table(sim_repertoire$j_call))

IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.016 0.032 0.074 0.346 0.122 0.410 

I tried to with number_of_seqs of 5000,10.000, 50.000 , but the results are the same.

Thanks in advance for the support :)

P.S: I also tried to change the freq_update_time paramter, but the result changes only slightly

freq_update_time 300
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.010 0.034 0.062 0.350 0.130 0.414

freq_update_time 350
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.024 0.040 0.068 0.326 0.146 0.396 

freq_update_time 400
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.020 0.036 0.060 0.350 0.142 0.392 

freq_update_time 200
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.016 0.028 0.068 0.328 0.140 0.420

freq_update_time 150
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.016 0.028 0.062 0.342 0.130 0.422 

freq_update_time 100
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.012 0.034 0.058 0.342 0.130 0.424

freq_update_time 50
IGHJ1 IGHJ2 IGHJ3 IGHJ4 IGHJ5 IGHJ6 
0.024 0.018 0.062 0.340 0.136 0.420 

Thanks for developing this amazing package to simulate repertoire sequences! I hope that you could provide a brief example of how to specify the the "chain" parameter to simulate paired chains?

Thanks,
Yuanyuan