text:
Working with Smart-seq3 data, have some questions regarding UMI and bar coding
Most of my experience has been working with genomic data, this time I am working on data from Smart-seq3. RNA excreted by cells in culture medium were processed with Smart-seq3 (From the looks of it, to capture as many fragments as possible). I have to perform differential expression analysis and check for any significant genes. I have the fastq files and trimmed them for adapters and polyA using trim_galore.
The protocol suggested processing fastq file using zUMIs, however the data isn't traditional single cell sequencing, so as per one comment of their procotol, they replied that for bulk sequencing-esque applications, we can just trim out the UMIs.
I have not found a decent way to remove UMIs from fastq file, are there any tools to remove them? From what I've read, UMI_tools can extract the UMI information from the fastq file, but how I remove it from the fastq file?

May I ask why Smart-seq3 protocol was used if you plan not to use UMIs? I understand you will want to trim away the cell barcodes (if the RT primer contained them), as they are meaningless in bulk sequencing, but UMIs are there to enable more accurate counting of RNA molecules, especially if the amount of starting material is low, which requires a high number of PCR cycles in library preparation.

More like, I have entered the project during the analysis stage, so I was not involved in designing the project, and now I have been tasked to analyze the data in the most meaningful way I personally can. edit2add - I have not worked with UMIs or RNA data. Only 30-40% of reads have UMI, so even after collapsing reads normalization and analysis might be an issue. Should I analyze UMI reads separately?

I have never worked with Smart-seq data, but in our lab we use a Drop-seq-like protocol for multiplexing bulk RNA-seq samples, and we process the data using the Drop-seq-tools software bundle. One idea could be to just use zUMIs as recommended for single-cell data, since it can be set up to handle Smart-seq v3 reads. Once you have the differential gene expression matrix, you can then continue with regular RNA-seq analysis tools.

Removing the UMIs may be trickier because, as you say, not all reads contain them. It looks like zUMIs can be set to output demultiplexed BAM files per cell (in your case per sample), so if you don't want to use UMIs, you could count the reads from those files, ignoring the UMI tags.


text:
Old program make error ProtDeform
Im trying to run ProtDeform on ubuntu 20.04
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3125161/
I need scoring function :/ to my project
But error occurs:
(base) marcin@marcin-B550-AORUS-ELITE-AX-V2:~/ProtDeformV1.0$ ./install
rm: cannot remove 'ProtDeform': No such file or directory
Compiling...
gcc -c  -g Matras.c Amino.c Vector.c GotohDPdis.c DiffStr.c BetaTopo.c RegressL.c Jacobi3.c SecStrEle.c ScoreENV.c ScoreDIS.c ScoreSSE.c BuildUpL.c DiheDis.c RMS.c OutAln.c StrLibL.c OutLibL.c LlistIO.c CalZsc.c LlistAlign.c MergeSort.c PairAli.c IO.c
Matras.c:77:1: warning: return type defaults to ‘int’ [-Wimplicit-int]
   77 | main(argc,argv)
      | ^~~~
Matras.c: In function ‘main’:
Matras.c:126:3: warning: implicit declaration of function ‘Set_AAnum’ [-Wimplicit-function-declaration]
  126 |   Set_AAnum();
      |   ^~~~~~~~~
Matras.c:390:27: warning: zero-length gnu_printf format string [-Wformat-zero-length]
  390 |         sprintf(envscfile,"");
      |                           ^~
Matras.c:431:23: warning: zero-length gnu_printf format string [-Wformat-zero-length]
  431 |   sprintf(PAR.COMMAND,"");
      |                       ^~
Matras.c:436:29: warning: implicit declaration of function ‘Get_Date_String’; did you mean ‘Get_Time_String’? [-Wimplicit-function-declaration]
  436 |   printf("#<DATE> \"%s\"\n",Get_Date_String());
      |                             ^~~~~~~~~~~~~~~
      |                             Get_Time_String
Matras.c:436:22: warning: format ‘%s’ expects argument of type ‘char *’, but argument 2 has type ‘int’ [-Wformat=]
  436 |   printf("#<DATE> \"%s\"\n",Get_Date_String());
      |                     ~^      ~~~~~~~~~~~~~~~~~
      |                      |      |
      |                      char * int
      |                     %d
Matras.c:442:10: warning: too many arguments for format [-Wformat-extra-args]
  442 |   printf("# QueFile \"%s\" LibFile \"%s\"\n",MODE,QueFile,LibFile);
      |          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Matras.c:498:4: warning: implicit declaration of function ‘Pairwise_3D_Alignment’ [-Wimplicit-function-declaration]
  498 |    Pairwise_3D_Alignment(&proA,&proB,&DmapA,&DmapB,&AHead,&SHead,inalifile,inssefile);
      |    ^~~~~~~~~~~~~~~~~~~~~
Matras.c:502:4: warning: implicit declaration of function ‘mainPD’; did you mean ‘main’? [-Wimplicit-function-declaration]
  502 |    mainPD(&proA, &proB, AHead.next, inalifile ); /* JR ProtDeform */
      |    ^~~~~~
      |    main
Matras.c:529:5: warning: implicit declaration of function ‘Output_Path_Horizontal_ClustalW’; did you mean ‘Output_Path_Horizontal_HTMLcolor’? [-Wimplicit-function-declaration]
  529 |     Output_Path_Horizontal_ClustalW(&proA,&proB,an,outclwfile);
      |     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
      |     Output_Path_Horizontal_HTMLcolor
Matras.c:726:6: warning: implicit declaration of function ‘Cal_Raw_MeanSD_Llist’; did you mean ‘Cal_Raw_Zscore_Llist’? [-Wimplicit-function-declaration]
  726 |      Cal_Raw_MeanSD_Llist(&ALibHead,&Mean,&SD);
      |      ^~~~~~~~~~~~~~~~~~~~
      |      Cal_Raw_Zscore_Llist
Matras.c:512:26: warning: ‘sprintf’ may write a terminating nul past the end of the destination [-Wformat-overflow=]
  512 |     sprintf(OutHead,"%s%d",OheadStr,an->num);
      |                          ^
Matras.c:512:5: note: ‘sprintf’ output between 2 and 139 bytes into a destination of size 128
  512 |     sprintf(OutHead,"%s%d",OheadStr,an->num);
      |     ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Matras.c:587:43: warning: ‘%s’ directive writing up to 127 bytes into a region of size between 113 and 240 [-Wformat-overflow=]
  587 |      { sprintf(command,"perl %s/rotSTR.pl %s.mat %s -nc A -S B > %s.pdb",
      |                                           ^~
  588 |              BASE_DIR,OutHead,proAfile,OutHead);
      |                       ~~~~~~~
Matras.c:587:8: note: ‘sprintf’ output between 40 and 484 bytes into a destination of size 256
  587 |      { sprintf(command,"perl %s/rotSTR.pl %s.mat %s -nc A -S B > %s.pdb",
      |        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  588 |              BASE_DIR,OutHead,proAfile,OutHead);
      |              ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Matras.c:590:62: warning: ‘%s’ directive writing up to 127 bytes into a region of size between 33 and 223 [-Wformat-overflow=]
  590 |        sprintf(command,"perl %s/rotSTR.pl - %s -nc B -S B >> %s.pdb",BASE_DIR,proBfile,OutHead);
      |                                                              ^~                        ~~~~~~~
Matras.c:590:8: note: ‘sprintf’ output between 38 and 355 bytes into a destination of size 256
  590 |        sprintf(command,"perl %s/rotSTR.pl - %s -nc B -S B >> %s.pdb",BASE_DIR,proBfile,OutHead);
      |        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Matras.c:594:43: warning: ‘%s’ directive writing up to 127 bytes into a region of size between 113 and 240 [-Wformat-overflow=]
  594 |      { sprintf(command,"perl %s/rotSTR.pl %s.mat %s -nc A -S P -ch %c -H W > %s.pdb",
      |                                           ^~
  595 |            BASE_DIR,OutHead,proAfile,proA.chain,OutHead);
      |                     ~~~~~~~
Matras.c:594:8: note: ‘sprintf’ output between 51 and 495 bytes into a destination of size 256
  594 |      { sprintf(command,"perl %s/rotSTR.pl %s.mat %s -nc A -S P -ch %c -H W > %s.pdb",
      |        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  595 |            BASE_DIR,OutHead,proAfile,proA.chain,OutHead);
      |            ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Matras.c:597:73: warning: ‘%s’ directive writing up to 127 bytes into a region of size between 23 and 213 [-Wformat-overflow=]
  597 |       sprintf(command,"perl %s/rotSTR.pl - %s -nc B -S P -ch %c -H W>> %s.pdb",
      |                                                                        ^~

  598 |         BASE_DIR,proBfile,proB.chain,OutHead);
      |                                      ~~~~~~~

Matras.c:597:8: note: ‘sprintf’ output between 48 and 365 bytes into a destination of size 256
  597 |        sprintf(command,"perl %s/rotSTR.pl - %s -nc B -S P -ch %c -H W>> %s.pdb",
      |        ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
  598 |          BASE_DIR,proBfile,proB.chain,OutHead);
      |          ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Amino.c: In function ‘Amino_Name’:
Amino.c:104:11: warning: implicit declaration of function ‘strncmp’ [-Wimplicit-function-declaration]
  104 |       if (strncmp(t,"ALA",3)==0)  r = 'A';
      |           ^~~~~~~
Amino.c:104:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  104 |       if (strncmp(t,"ALA",3)==0)  r = 'A';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:105:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  105 |  else if (strncmp(t,"GLU",3)==0)  r = 'E';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:106:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  106 |  else if (strncmp(t,"GLN",3)==0)  r = 'Q';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:107:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  107 |  else if (strncmp(t,"ASP",3)==0)  r = 'D';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:108:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  108 |  else if (strncmp(t,"ASN",3)==0)  r = 'N';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:109:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  109 |  else if (strncmp(t,"LEU",3)==0)  r = 'L';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:110:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  110 |  else if (strncmp(t,"GLY",3)==0)  r = 'G';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:111:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  111 |  else if (strncmp(t,"LYS",3)==0)  r = 'K';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:112:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  112 |  else if (strncmp(t,"SER",3)==0)  r = 'S';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:113:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  113 |  else if (strncmp(t,"VAL",3)==0)  r = 'V';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:114:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  114 |  else if (strncmp(t,"ARG",3)==0)  r = 'R';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:115:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  115 |  else if (strncmp(t,"THR",3)==0)  r = 'T';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:116:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  116 |  else if (strncmp(t,"PRO",3)==0)  r = 'P';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:117:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  117 |  else if (strncmp(t,"ILE",3)==0)  r = 'I';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:118:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  118 |  else if (strncmp(t,"MET",3)==0)  r = 'M';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:119:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  119 |  else if (strncmp(t,"PHE",3)==0)  r = 'F';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:120:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  120 |  else if (strncmp(t,"TYR",3)==0)  r = 'Y';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:121:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  121 |  else if (strncmp(t,"CYS",3)==0)  r = 'C';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:122:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  122 |  else if (strncmp(t,"TRP",3)==0)  r = 'W';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:123:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  123 |  else if (strncmp(t,"HIS",3)==0)  r = 'H';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:124:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  124 |  else if (strncmp(t,"UNK",3)==0)  r = 'u';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:125:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  125 |  else if (strncmp(t,"ACE",3)==0)  r = 'a';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
Amino.c:126:27: warning: ‘strncmp’ argument 3 type is ‘int’ where ‘long unsigned int’ is expected in a call to built-in function declared without prototype [-Wbuiltin-declaration-mismatch]
  126 |  else if (strncmp(t,"FOR",3)==0)  r = 'f';
      |                           ^
<built-in>: note: built-in ‘strncmp’ declared here
GotohDPdis.c: In function ‘GotohDP_Dscore_Iterative’:
GotohDPdis.c:202:3: warning: implicit declaration of function ‘ScoreDistanceSPair’; did you mean ‘ScoreDistanceS’? [-Wimplicit-function-declaration]
  202 |   ScoreDistanceSPair(proA,proB,an,DmapA,DmapB);
      |   ^~~~~~~~~~~~~~~~~~
      |   ScoreDistanceS
GotohDPdis.c: In function ‘GotohDP_Dscore’:
GotohDPdis.c:223:33: warning: format ‘%s’ expects argument of type ‘char *’, but argument 2 has type ‘int’ [-Wformat=]
  223 |   { printf("#ERROR:PRE.MaxNali %s is less than Nali %d\n",PRE.MaxNali,an->Nali);
      |                                ~^                         ~~~~~~~~~~~
      |                                 |                            |
      |                                 char *                       int
      |                                %d
DiffStr.c: In function ‘Number_Sec’:
DiffStr.c:202:5: warning: type of ‘sym’ defaults to ‘int’ [-Wimplicit-int]
  202 | int Number_Sec(sym,Nstr)
      |     ^~~~~~~~~~
SecStrEle.c: In function ‘Output_SSE_PDB’:
SecStrEle.c:355:13: warning: too many arguments for format [-Wformat-extra-args]
  355 |  fprintf(fp,"REMARK tFactor : 0 for Cend of SSE, 200 for  helix, 100 for sheet\n",pr->pname);
      |             ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
ScoreENV.c: In function ‘Number_Sec’:
ScoreENV.c:313:5: warning: type of ‘sym’ defaults to ‘int’ [-Wimplicit-int]
  313 | int Number_Sec(sym,Nstr)
      |     ^~~~~~~~~~
ScoreENV.c: In function ‘Read_Env_Score’:
ScoreENV.c:264:22: warning: ‘sprintf’ may write a terminating nul past the end of the destination [-Wformat-overflow=]
  264 |  sprintf(Fname,"%s/%s",SCORE_DIR,fname);
      |                      ^
ScoreENV.c:264:2: note: ‘sprintf’ output 2 or more bytes (assuming 129) into a destination of size 128
  264 |  sprintf(Fname,"%s/%s",SCORE_DIR,fname);
      |  ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BuildUpL.c: In function ‘Build_Up_LL’:
BuildUpL.c:256:2: warning: implicit declaration of function ‘Merge_Sort_Double_Linked_List’ [-Wimplicit-function-declaration]
  256 |  Merge_Sort_Double_Linked_List(head,'R');
      |  ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
BuildUpL.c: In function ‘Read_SSE_Alignment’:
BuildUpL.c:677:4: warning: implicit declaration of function ‘Get_Part_Of_Line’ [-Wimplicit-function-declaration]
  677 |    Get_Part_Of_Line(num1,word1,1,strlen(word1));
      |    ^~~~~~~~~~~~~~~~
OutAln.c: In function ‘Write_Rasmol_Color_Alignment’:
OutAln.c:564:34: warning: implicit declaration of function ‘Get_Date_String’ [-Wimplicit-function-declaration]
  564 |  fprintf(fp,"# DATE      '%s'\n",Get_Date_String());
      |                                  ^~~~~~~~~~~~~~~
OutAln.c:564:28: warning: format ‘%s’ expects argument of type ‘char *’, but argument 3 has type ‘int’ [-Wformat=]
  564 |  fprintf(fp,"# DATE      '%s'\n",Get_Date_String());
      |                           ~^     ~~~~~~~~~~~~~~~~~
      |                            |     |
      |                            |     int
      |                            char *
      |                           %d
OutAln.c: In function ‘Write_Rasmol_Color_Alignment_Each_Protein’:
OutAln.c:715:27: warning: format ‘%s’ expects argument of type ‘char *’, but argument 3 has type ‘int’ [-Wformat=]
  715 |  fprintf(fp,"# DATE      %s\n",Get_Date_String());
      |                          ~^    ~~~~~~~~~~~~~~~~~
      |                           |    |
      |                           |    int
      |                           char *
      |                          %d
OutAln.c:731:10: warning: implicit declaration of function ‘Ncompare’ [-Wimplicit-function-declaration]
  731 |  Ncomp = Ncompare(a);
      |          ^~~~~~~~
OutAln.c: In function ‘Write_Molscript_Alignment_Each_Protein’:
OutAln.c:880:27: warning: format ‘%s’ expects argument of type ‘char *’, but argument 3 has type ‘int’ [-Wformat=]
  880 |  fprintf(fp,"! DATE      %s\n",Get_Date_String());
      |                          ~^    ~~~~~~~~~~~~~~~~~
      |                           |    |
      |                           |    int
      |                           char *
      |                          %d
OutAln.c:905:35: warning: implicit declaration of function ‘Remove_Space’ [-Wimplicit-function-declaration]
  905 |        if (OutProteinType == 'A') Remove_Space(rsta,pA->Rnum[a->nalA[i]]);
      |                                   ^~~~~~~~~~~~
OutAln.c: In function ‘Find_nal’:
OutAln.c:962:5: warning: type of ‘q’ defaults to ‘int’ [-Wimplicit-int]
  962 | int Find_nal(q,Nali,nalA)
      |     ^~~~~~~~
OutAln.c: In function ‘Output_Alignment_Information’:
OutAln.c:1058:2: warning: implicit declaration of function ‘Get_Part_Of_Line’ [-Wimplicit-function-declaration]
 1058 |  Get_Part_Of_Line(ln,pA->compnd,0,40);
      |  ^~~~~~~~~~~~~~~~
StrLibL.c: In function ‘Read_Library_List’:
StrLibL.c:75:5: warning: implicit declaration of function ‘Get_First_Field’ [-Wimplicit-function-declaration]
   75 |     Get_First_Field(name,tax,buffer);
      |     ^~~~~~~~~~~~~~~
StrLibL.c: In function ‘Check_Library_List’:
StrLibL.c:113:3: warning: implicit declaration of function ‘Get_Part_Of_Line’ [-Wimplicit-function-declaration]
  113 |   Get_Part_Of_Line(dir,bn->Name,1,2);
      |   ^~~~~~~~~~~~~~~~
OutLibL.c: In function ‘Output_Library_ClustalW’:
OutLibL.c:234:4: warning: implicit declaration of function ‘Get_Part_Of_Line’ [-Wimplicit-function-declaration]
  234 |    Get_Part_Of_Line(ln,bn->sim->Name,0,16);
      |    ^~~~~~~~~~~~~~~~
LlistIO.c: In function ‘Read_BSSP_File_LinearList’:
LlistIO.c:161:2: warning: implicit declaration of function ‘Set_SSE’ [-Wimplicit-function-declaration]
  161 |  Set_SSE(pr);
      |  ^~~~~~~
LlistIO.c:162:2: warning: implicit declaration of function ‘Set_Dis_bwn_SSE’ [-Wimplicit-function-declaration]
  162 |  Set_Dis_bwn_SSE(pr,'C');
      |  ^~~~~~~~~~~~~~~
LlistIO.c: In function ‘Translate_BSSP_File_LinearList’:
LlistIO.c:218:51: warning: implicit declaration of function ‘tolower’ [-Wimplicit-function-declaration]
  218 |             if (X->line[75]=='C') pr->Seq[Rnum] = tolower(pr->Seq[Rnum]);  /* for CIS conf */
      |                                                   ^~~~~~~
PairAli.c: In function ‘Pairwise_3D_Alignment’:
PairAli.c:188:3: warning: implicit declaration of function ‘Set_Value_AlignNode’; did you mean ‘Tail_Value_AlignNode’? [-Wimplicit-function-declaration]
  188 |   Set_Value_AlignNode(Ahead,'S');    /* JR  It was  D */
      |   ^~~~~~~~~~~~~~~~~~~
      |   Tail_Value_AlignNode
PairAli.c:189:3: warning: implicit declaration of function ‘Merge_Sort_Double_Linked_List’ [-Wimplicit-function-declaration]
  189 |   Merge_Sort_Double_Linked_List(Ahead,'R');
      |   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~
IO.c: In function ‘Read_BSSP_File’:
IO.c:705:3: warning: implicit declaration of function ‘Set_SSE’ [-Wimplicit-function-declaration]
  705 |   Set_SSE(pr);
      |   ^~~~~~~
IO.c:706:3: warning: implicit declaration of function ‘Set_Dis_bwn_SSE’ [-Wimplicit-function-declaration]
  706 |   Set_Dis_bwn_SSE(pr,'C');
      |   ^~~~~~~~~~~~~~~
IO.c: In function ‘Read_Sequence’:
IO.c:886:13: warning: implicit declaration of function ‘isalpha’ [-Wimplicit-function-declaration]
  886 |    else if (isalpha(sym)!=0)
      |             ^~~~~~~
IO.c: In function ‘Read_Matras_Environment_File’:
IO.c:87:23: warning: ‘%s’ directive writing up to 63 bytes into a region of size between 0 and 127 [-Wformat-overflow=]
   87 |   sprintf(Matfile,"%s/%s",homedir,matfile);
      |                       ^~          ~~~~~~~
IO.c:87:3: note: ‘sprintf’ output between 2 and 192 bytes into a destination of size 128
   87 |   sprintf(Matfile,"%s/%s",homedir,matfile);
      |   ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
IO.c: In function ‘Read_PDB_File’:
IO.c:378:20: warning: ‘%s’ directive writing up to 3 bytes into a region of size between 0 and 127 [-Wformat-overflow=]
  378 |  sprintf(Fname,"%s/%s/%s",PDB_DIR,dir,fname);
      |                    ^~             ~~~
IO.c:378:2: note: ‘sprintf’ output 3 or more bytes (assuming 133) into a destination of size 128
  378 |  sprintf(Fname,"%s/%s/%s",PDB_DIR,dir,fname);
      |  ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
g++ -c  -g vector.cpp TransSeq.cpp Geometry.cpp WriteDeformPDB.cpp Evaluation.cpp  -w
Geometry.cpp:30:31: error: ‘<unnamed enum> GEOMETRY’, declared using unnamed type, is used but never defined [-fpermissive]
   30 | extern enum { RIGID, AFFINE } GEOMETRY;
      |                               ^~~~~~~~
WriteDeformPDB.cpp:16:31: error: ‘<unnamed enum> GEOMETRY’, declared using unnamed type, is used but never defined [-fpermissive]
   16 | extern enum { RIGID, AFFINE } GEOMETRY;
      |                               ^~~~~~~~
make: *** [Makefile:125: ../ProtDeform] Error 1

I've contacted the program developers about this, but am still waiting for answer from them.
Can you help?
EDIT1:
ox
I've checked
Ubuntu 16.04 on virtualbox
Ubuntu 20.04 with i7
Ubuntu 20.04 with Ryzen 9
same error...
EDIT 2:
I have found hidden file .matras:
###############################
### MATRAS ENVIRONMENT FILE ###
###############################

BASE_DIR .
SCORE_DIR  data_sc/ROM-04JAN29
TMP_DIR    .
BSSP_DIR  /DB/BSSP
PDB_DIR   /DB/PDB

I'm checking to see if it simply requires PDB

What is ProtDeform? Why are you trying to use it? Can you give any additional context around your problem that would help in finding a solution?

@gringer I add some information. Yes, still waiting for answer :/

The tool is more than 10 years old, without documentation for compilation and without any documentation towards software versions. The paper has a single citation, indication it never has been adopted by the community. Hence, find something else that actually works. Not what you want to hear, I know.

is this webserver giving any result with your input https://bioinfo.uib.es/~recerca/ProtDeform/

are you trying to run or install the software ??

@pippo1980 install to run :)

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
paired-end short reads: will one file suffice?
I have a quick question about paired-end short reads. I have multiple genomes that were sequenced with paired-end Illumina NextSeq 200 technology, resulting in two fastq files per sample: 001_R1.fastq.gz and 001_R2.fastq.gz, typical for short read data.
I would like to use a k-mer counting and clustering approach to put these genomes together on a phylogenetic tree. I was able to use both reads in one go to count k-mers but, as you can probably guess, this resulted in a very high frequency of k-mers identified and accounted for.
For this specific use case (not interested in changing my approach), would it make more sense to just use a single file for kmer counting (i.e. only the 001_R1.fastq.gz file)?
EDIT:
For more context, I'm trying to do an alignment-free approach to produce a phylogenetic distance comparison of these strains. I was using long reads before and it worked fine but Nanopore reads tend to be very messy. I have short-read data for all of the strains and would like to use that instead but that comes at the price of an overload of data.
I suppose my main concern/question is how do I pre-process these short reads without assembling or aligning to a reference to count kmers? Or is that not possible without at least alignment?

Do you mean NextSeq 2000?

I would be concerned counting them twice. What you could do is map the paired end reads then count kmers. Or you could stitch the overlapping R1 and R2. Or just do not change what you are doing. If accurate mapping is important, I would suggest to map the reads base don R1 and R2 then count kmers without duplicating the count.

Each of the regions sequenced will be covered by tens or hundreds or even thousands of reads depending on coverage. How are you taking that into account? If you are treating the reads as independent sequences, you will be over counting by orders of magnitude whether you use one file or both.

Would an alignment be the only way to preprocess these short reads for analysis? Ideally, it would be a completely alignment-free approach but do you mean align to a reference or align the R1.fastq to the R2.fastq? @Supertech

That's what I'm asking here @terdon ... I honestly don't know how to process short reads for this kind of approach. I know there are a lot of reads here because it's very high coverage but I also know I have been asked to not use an assembly because we want to conserve haplotypes and to not use an alignment to a reference to not bias the seqeunce...

I suspect you won't be able to avoid assembling or aligning. If you are working with strains that don't have their own, well studied reference genome and so want to avoid aligning so as not to bias your data towards whichever strain has been chosen as the reference for your species (please correct me if this isn't what's going on), then you will need to assemble a denovo genome from these reads. At the very least you will need to build the longest consensus contigs you can. Unfortunately, this isn't something I have experience with so I can't help you do this beyond telling you that, as I'm sure you already know, there are tools out there that will try to assemble your reads into contigs.
I will instead try to explain why I think you can't avoid assembly. Short read fastq files, paired or not, contain short sequences of ~150-200 base pairs each. Now, depending on the depth of coverage, each region of your genome will be covered by multiple such short reads. Therefore, the data you have in your fastq files will be extremely repetitive and redundant since each of your genome's regions will be represented by multiple reads, and any attempt to find k-mers will be dogged by this redundancy. This is why you got so many more results than you expected and I don't see how using one file only would help at all.
If you had sequenced using an assay that includes unique molecular identifiers (UMIs), then you would be able to extract non-redundant, representative reads for all regions but I cannot think of a way to do this without UMIs off the top of my head. Not saying it isn't possible, only that you need someone more clever than I.
That said, the best I can come up with is that if your samples were sequenced in the same sequencing run and so can be assumed to have the same sequencing biases, you might be able to build your trees despite the redundancy in the data. My thinking here is that what you are really interested in isn't the k-mers themselves or their level of over-representation in each species, but only the difference in the k-mer content across your samples. If so, given that you sequenced all samples in the same run, they should all have the same level of over-representation and any differences might be enough for you to be able to cluster/build your trees.
This seems worth a try to me, but I haven't worked on something like this in many years so if you go down this route, make sure to double check and ensure I'm not leading you astray!

The problem with using single-end short reads is that they're often not long enough to uniquely match to a database (especially for highly-conserved genes, like the ribosomal RNAs). Whether or not this is a problem for your particular situation is unclear; you haven't provided enough context / information for me to work that out. For example, a single end of a 2x50bp run on a NextSeq 2000 might cause problems (especially for medium to large eukaryote genomes), but a 2x300bp run on a NextSeq 2000 is probably okay.
FWIW, Kraken2 has a paired mode, which combines the results from both read ends from the same template source:

Usage of --paired also affects the --classified-out and --unclassified-out options; users should provide a # character in the filenames provided to those options, which will be replaced by kraken2 with "_1" and "_2" with mates spread across the two files appropriately. For example:

kraken2 --paired --classified-out cseqs#.fq seqs_1.fq seqs_2.fq

See the Wiki / manual here:
https://github.com/DerrickWood/kraken2/wiki/Manual#classification

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Predict population based on PC coordinates
Utilizing a reference panel, I want to assign most likely population label to each individual in the study. Following are the files I have:

Reference panel population labels:
HG00121 EUR
HG00122 EUR
HG00123 EUR

Reference panel PC (Principal Component) coordinates:
HG00121 -0.0108487      0.0259052       -0.0099605      0.0168545
HG00122 -0.0115789      0.0264642       -0.01722        0.0157885
HG00123 -0.00996687     0.0266983       -0.0155755      0.0155963

Study individual PC coordinates:
Study_individual1 -0.0108487      0.0259052       -0.0099605      0.0168545

How do I go about assigning a population label to the study individual based of the reference panel? Most of the softwares I have come across such as evaluate_check_ancestry(plinkQC), FRAPOSA rely on inputs provided to PLINK which is not helpful in my case.

I meant Principal Components (PC) coordinates. Edited question.

Why start with PCA? It's a means of reducing dimensionality to retain the highest amount of variability, but those reduced dimensions may have nothing at all to do with differences between the categories. Suppose you find that the query study is very similar to Study X in PC1 and PC2 - that doesn't necessarily imply it's the same category as Study X, as PC1 and PC2 may totally irrelevant to the category.

@NuclearHoagie In the era of machine learning PCA is a common starting point. It's classed as unsupervised learning and used to generate the hypothesis for supervised learning. Also in the world of RNAseq (which this isn't) its frontline approach. In pop gen ... generally its one of a number of starting points, but its not the end analysis.

@M__ Right, I just don't think it's a particularly good starting point for this problem. If you want to predict some label based on feature values, you probably shouldn't start by throwing away most of the data without having any idea if what you're discarding is useful for prediction or not. If study label is easily determined by a low-variance feature surrounded by lots of high-variance noise, using only the top PCs will make the task all but impossible, when it could be easily accomplished by starting with a supervised method. There's no reason why "distance" in PC space should be "better".

@NuclearHoagie thats an interesting way of expressing a dimensionality reduction critique. I agree on the heterogeneity issue - it's definitely an issue and weakness - depending how the variance is spread across PCs, and needs further analysis. I'll be honest I start with supervised methods too, but in any walk of life someone often deploys unsupervised analysis.

@NuclearHoagie it turns out that for human genetic data, the major principals components do align with major population groups and lower PCs strongly correlate with finer scale populations - it’s not an ‘impossible task’ at all. PCA is fine for this purpose if you need to do it quickly for lots of samples, particularly if you are just interested in continental-scale groups like the OP seems to be interested in

I apologise for being negative about this, but I am afraid you cannot use principle components for population inference other than assigning clusters within a bivariate or 3D graph. Eigen values are complex derivatives of data to achieve dimensionality reduction.
For population genetics you require a sample size => 10. I don't know the organism nor the question but goal is to derive allelic frequencies for each genome position. A common metric is F-statistics, for example Fst or Hardy-Weinberg equilibrium.

I do have study sample size of 10 (tried to make the question by including just 1). Thank you for the input. I will look into Fst and Hardy-Weinberg equilibrium

Dear @user2998764 upvotes/accepts much preferred to thanks, but thanks. I'm glad this helped.  10 is an okay sample size and you can always ask further given your experimental aims for the most useful metric.

‘You cannot use PCA for population inferrence other than assigning clusters’ but isn’t that basically what the OP is trying to do?

Its application is for the data mining phase, not the final answer. At a minimum the OP would need to assess cumulative variance across PCs and it gets complicated.

Here is a very simple way to get your started. Of course you can make it more complex:
Take the mean PC1 and PC2 coordinates for each 'reference panel population label' and then calculate the reference panel population which has the closest Cartesian distance to your target, i.e. $\sqrt{{(PC1_{r} - PC1_{t})^2 + (PC2_{r} - PC2_{t})^2}}$, where $PC1_{r}$ is the mean PC1 for a ref population and $PC1_{t}$ is PC1 for the target.


text:
Identifying the promoter on the reverse strand using visualization software (seqmonk)
I feel like this might be a really silly question and I've been going back and forth for awhile and I think confusing myself- so sorry if this is dumb.
I am using Seqmonk to visualize genes in an assembly (but I suspect the issue would be the same in any visualization software). I am trying to identify promoters of specific genes to determine if they are differentially methylated. Genes on the forward strand (highlighted in red), the promoter would inherently be upstream- in this case to the left (i.e. if my start site was 5000 and the end site 6000, I would assume I could find my promoter to the left of the start site between 3000 and 5000 bases.
My confusion is what happens with genes that are on the reverse strand (in blue here). Start site and end site could be the same (5000 and 6000- as end sites are always bigger than the start site). Does this mean that in the software, the promoter is expected to be on the right hand side (i.e. bases 6000-8000)? This has been my going thought but now I wonder if having bases that seem downstream (because they're bigger) makes sense, even though it really does because downstream on the forward strand would be upstream on the reverse strand. Or if I need to always do 3000-5000 bases as the promoter, regardless of the the strand of origin.

If I'm looking at a gene on the reverse strand, where would I find the promoter.

If you intended to post an image, we are not able to see it. Not clear what red and blue are.

Have a look at this UCSC browser shot.

Looking at the WDR92 transcript, note that it is transcribed right to left. This would indicate that the promoter is on the right. Indeed, we see a variety of promoter histone marks on the right of the gene, but not on the left.
In other words, yes, for genes on the reverse strand, the larger coordinate will be at the beginning of the gene, and the promoter will include a larger coordinate than the rightmost transcribed base.

THANK YOU. This is what I thought but I suddenly got wrapped up in my head. Super helpful.

@user18170 No problem. consider upvoting/accepting the answer if you found it helpful, in case others find this in the future.


text:
How to get cytoband and gene level copy number from genome wide SNP array copy number data?
I have (human) Illumina genome wide SNP array copy number data. For each SNP genomewide, I have Log R Ratio (LRR) and B Allele Frequency (BAF).
What tool(s) can I use to get the integer copy numbers (either -3 to +3 or 0 to inf) for each cytoband and for each gene in the genome?
I assume that people routinely do this? Are there any R packages/pipelines/code to do this?
Do people use DNAcopy and/or GISTIC for this?
Thank you.

To obtain cytoband and gene-level copy numbers from genome-wide SNP array copy number data, you might want to consider using the ASCAT tool, as it makes allele-specific copy number segment data and derived integer gene-level copy number values from SNP array data. The SNP6 copy number analysis pipeline (ASCAT2) I believe is specifically for this purpose.
See also:

CNVkit, which has a method for copy number detection via using targeted reads and the non-specifically captured off-target reads to infer copy number evenly across the genome.

Array Comparative Genomic Hybridization (aCGH) and SNP Genotyping Technologies

CNVannotator

Thank you very much for your answer with multiple options. I will look into these. ¡Gracias! :)


text:
I am trying to create a subset of 10k variants from 25-30 unmapped contigs of a g.vcf file including the header
My objective is to take a g.vcf.gz file and from 25-30 unmapped contigs with titles like "NW_020192317.1", I want to make a subset of ~10k variants from each of the unmapped contigs and make one final g.vcf file that includes the header from the original g.vcf.gz file.
From the post here, I have tried running:
samtools view -bo output.g.vcf -s 123.4 inputfile.g.vcf.gz NW_020192317.1

and I keep receiving the output:
Aborted (core dumped)

I'm not too picky about the number of variants per contig, as long as it is in the thousands, that's why I tried using ".4"
Alternatively, I tried using the code:
samtools view input.g.vcf.gz "NW_020192317.1:1-2" > output.g.vcf

and I recieved the same response of "Aborted (core dumped)" with nothing else. I am confused what is wrong here and why it is throwing an error.

You at least need to be using bcftools, not samtools, if you are working with vcf files. samtools is designed for bam/sam, not vcf, and that’s why it’s throwing that error

Do you have a list of contigs of interest? Are those contig names the first field of the VCF? Does this need to be a random selection or can you keep the first 10k variants from each contig?

Also, gvcf files contain both variants non-variant lines. Do you only want variants?

After some guidance from my advisor, this command allowed me to get the output file I needed.
bcftools view -o newvcfFileName.g.vcf -s 123.2 locationOfTheFile chrom1 chrom2 ... chrom23(mito DNA) --force-samples

I listed every chromosome's title until I typed --force-samples. I found it simpler to copy and paste the names in.


text:
How can I develop a snakemake rule for all files in multiple directories?
I am a beginner in snakemake and I'm trying to develop a working flow.
Imagine I have multiple folders and in each folder I have a tumour and normal sample files:
SRC120/
      SRC120-T1.fastq
      SRC120-N1.fastq

SRC121/
      SRC121-T1.fastq
      SRC121-N1.fastq

and so on (in total I have more than 100 SRC folders).
I am trying to make a snakemake rule to make a symlink of these fastq files.
What I have written so far is:
import os

SAMPLES=os.listdir("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/")

def getFastqFilesForTissu(wildcards):
    fastqs = list()
    # Beware no other files than fastqs should be there
    for s in os.listdir("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/"+wildcards.sample):
        fastqs.append(os.path.join("/analysis/Anna/snakemake-demo/data/samples_fastq-test1/",wildcards.sample,s))
    return fastqs

rule all:
    input:
        expand("results/{sample}/{sample}.fastq", sample=SAMPLES)

rule foo:
    input:
        fastq_files = getFastqFilesForTissu
    output:
        "results/{sample}/{sample}.fastq"
    shell:
        'ln -s {input} {output}'

This works if I have one file in each folder and when the file name matches with the folder's name (i.e., when there is no T1 or N1 suffix).
Is there any way I could modify my rule to get it working for all samples in SRC directories? Or a better way to handle this issue?

What is the output you expect? What values do you want your {input} and {output} to have for a given example? This looks like a simple shell operation, do you really need to use snakemake?

This is just a very small part of a big workflow. the output is symlinked fastq files in the results folder

What symlinked files? What name should the symlink have and what should its target be? The more detail you give, the better the chances someone will be able to answer. Please [edit] your question and add the expected result so we know what you need.

All of these information is in the code if you read through the snakemake code. There is no need for any edition.

I suggest using package glob, which supports unix-style pattern search. For instance finding all .py files in my repo:
import glob
all_python_files = glob.glob(r'c:\personal\coderepo\**\*.py', recursive=True)

Which outputs full paths:

'c:\personal\coderepo\SimpleHTR\data\checkDirs.py',
'c:\personal\coderepo\....

I typically like to input the glob.glob to Pandas dataframes:
all_files = glob.glob(f'{base_folder}/*')
all_files = pd.DataFrame(all_files, columns=['full_path'])
all_files['my_filename'] = all_files['full_path'].apply(lambda x: os.path.basename(x))

#do whatever filtering is required:
all_files = all_files[all_files['my_file_name'].str.contains('regex_here')]

Then perform what ever I need to do in loop:
for idx, _ in all_files.iterrows():
    do_something(_['full_path'])


text:
Genetic relationships between Cyanobacteria: terrestrial vs aquatic
I believe one of the site members is an expert in Cyanobacteria hence the question.
I'm looking for a comprehensive 'phylogenetic tree' of the phylum Cyanobacteria and want to identify all the terrestrial members/groups. For example, Nodosilinea is one of the terrestrial members.
Part of my rationale for the question is that the Orders can be polyphyletic. From an outsider's point of view this is a difficult taxonomic group to understand.
Orders

Chroococcales
Chroococcidiopsidales
Gloeobacterales
Nostocales
Oscillatoriales
Pleurocapsales
Spirulinales
Synechococcales
Incertae sedis

Background
Cyanobacteria is a massive phylum and important to both terrestrial and aquatic habitats, but what I need to understand is how terrestrial members are distributed across the tree. The reason they are interesting (both aquatic and terrestrial) is that they are efficient at removing CO2 from the atmosphere. Reducing a carbon footprint is an area of growing commercial interest.

this is Benjy. I am a grad student at University of Nevada, Las Vegas. I am interested in terrestrial cyanobacteria too. I have read all the papers that you have mentioned such as Chen et al., 2021 and Strunecky et al., 2023. I would like to ask if you have an account of Google Scholar or Twitter. I am also working in a paper related to terrestrial cyanobacteria and bioinformatics and I'd like to read your future paper. Best,
Benjy

Yes however, I work in industry - a Big Five Tech is the key client (climate health) - and my company (not the client) protects its IP by non-disclosure. @Laura does publish in this area however. I rarely lift my anonymity, except that I that live in Mexico and I'm an active Academic Editor.

Ah! Finally, a question I can answer.
There has been some pretty good work about this recently. In particular, I would reccomend this paper by Chen et al. , which groups cyanobacteria into terrestrial / freshwater / marine and...well, the title says the rest. Some are edge cases by necessity (e.g. some coastal biofilm-forming cyanobacteria may be more similar to terrestrial ones than to many marine ones) but that's just nature.
Chen, MY., Teng, WK., Zhao, L. et al. Comparative genomics reveals insights into cyanobacterial evolution and habitat adaptation. ISME J 15, 211–227 (2021). https://doi.org/10.1038/s41396-020-00775-z (https://rdcu.be/dmJpm)
I regularly use this in my work, mainly because they have their trees posted on iTOL (in the 'Data availability' section) and this makes it really easy to search for specific strains.
Which brings us the the main problem. Right now there are about three main ways of separating the cyanobacterial phylum into orders. There is the NCBI database classification (largely based on 16S rRNA), like for example the Chen et al. paper uses. This may sometimes appear polyphyletic when overimposed on a genome tree. Then there is the GTDB classification (https://gtdb.ecogenomic.org/), which is based on genomes. This can be seen in Alcorta et al. 2020 (Alcorta J, Alarcón-Schumacher T, Salgado O and Díez B (2020) Taxonomic Novelty and Distinctive Genomic Features of Hot Spring Cyanobacteria. Front. Genet. 11:568223. doi: 10.3389/fgene.2020.568223) - on average, the tree shape is similar to the one in Chen et al., but the orders make more sense. The GTDB classification is good phylogenetically, and well-integrated with the database (for example, you can just get the list of all genomes/metagenomes available for an order). But because it's not focused on the taxonomy (the naming of things), the actual names may not always fit with the Bacterial/Botanical Code. One research group has been kind enough to go through the absolute mess that is cyanobacteria taxonomy, and now this paper exists.
Strunecký, O., Ivanova, A. P., & Mareš, J. (2023). An updated classification of cyanobacterial orders and families based on phylogenomic and polyphasic analysis. Journal of Phycology, 59(1), 12-51. (https://onlinelibrary.wiley.com/share/A5UXIHBUCIH9JRTIWFSE?target=10.1111/jpy.13304)
As far as I've seen, their orders more or less match with GTDB, just under different names (except for GTDB Cyanobacteriales with is split into more orders in this paper). They also have a much larger 16S rrNA phylogeny constrained by the genome one.
Hower, that solves only part of the problem. Because at the smaller scale (genus, species), you still have strains labeled by their morphologies, and culture collections are (for understandable reasons) hesitant to change their naming system. For particularly visually non-descript groups like 'Synechococcus' or 'Leptolyngbya', there is a good chance that a particular strain could be nearly anywhere in the phylogenetic tree. I would always get a 16S rRNA (if not already avaiable), BLAST it, and then check where related strains appear in a genome tree.
Edit: I have also heard of Cydrasil (https://www.cydrasil.org/), an online tool which can place a 16S rrNA sequence within the context of the cyanobacterial tree, but I haven't tested it yet.
Hope that helps!

A HUGE thanks. It's really helpful. It will take me a while to read through all this, but I can now have some expertise in cyanobacteria (or at least half competent). I work in industry BTW  (carbon neutral status for blue-chip clients) ... the analytics (and findings) are very advanced ... but my company will not publish (no comment)  ... My thinking is that cyanobacteria are huge players in the carbon cycle, but you'll know much more than I.

@M__ Hello.
Sorry to necro your question, just wanted to ask you something.
My work for the last few years is finally getting close to a publishable format, and I just wanted to ask whether you would agree to being acknowledged in my next paper. You've been answering my questions on Bioinformatics.SE for, well, years, and as someone who got into bioinformatics more deeply than my immediate group could help with, that was very valuable to me. I am very grateful, and would be glad to express that.

Hi @Laura, thats  cool and appreciated, thanks.


text:
Why does DNABert use overlapping k-mers as input?
I was reading through the DNABERT paper and found that their input features were k-mers. This is essentially using rolling window features.
Why train a transformer-like model on sequences of the form  ${a_n,..,a_{n+k}}$ rather than training it on the complete sequence ${a_1,...,a_{m}}$ with positional encoding?
To directly quote the paper:

Instead of regarding each base as a single token, we tokenized a DNA
sequence with the k-mer representation, an approach that has been
widely used in analyzing DNA sequences. The k-mer representation
incorporates richer contextual information for each deoxynucleotide
base by concatenating it with its following ones. The concatenation of
them is called a k-mer.

There doesn't seem to be any ablation study in this paper showing that this really is of any benefit.
I would appreciate any more information: I just don't see how the 3-mers  {ATG, TGG, GGC, GCT} contain more information than the sequence {A,T,G,G,C,T}, 1-mers, with positional encoding.

@M_ While DNABERT is in the title and is the only example used, it's only because I am quite new to bioinformatics. Moreover, the type of tokenization is completely irrelevant to my question. The question in the body is my main one; Why would you train any model with sequence of k-mers at all? Rolling window features makes sense in many contexts, but I cannot find any example anywhere of rolling window features being used the way they are in the DNABERT paper.

@M_  So, you're saying that k-mer for k>1 allows that but not for k=1? Please point me to the source saying that rolling windows are required for transformers to approximate non-linear functions. Which version of the universal approximation theorem is this?

@M_ Look, saying "this is really complicated" is not answering the question. Even listing a set of books I'll have to read to understand your claim would be better. You don't need k-mers to use language model; If you tokenize 1-mers, you essentially have a language with four words and can thus use language models. I've done this in proteomics before (having a 21 word language, one for each amino acid, in that case).

Let us continue this discussion in chat.

My answer below is about k-mers with k>1. For positional embedding, if you look at Fig. 1b you can see DNABERT did use positional encodings.

@Alexlok Actually, I just discovered there's a DNABERT-2 by the same team where they abandon k-mers altogether and get better performance (writing " However, we argue that the
computation and sample inefficiencies introduced by k-mer tokenization are primary obstacles in
developing large genome foundational models." in the abstract). So, I no longer believe this needs any explanation; Reading it, it feels like using k-mers was a mistake. Fair point about figure 1b, though; I completely missed that.

NLP and vectorisation is used to trap non-linearity which is highly applicable to sentences written in English.
Vectorisation of k-mers for an overlapping sliding window is weird because it's trying to assess the functionality of non-coding sequence, i.e. "silent" and non-cds mutations will be the dominant mutation class, but it would identify localised co-mutations, i.e. within a k-mer.
The only application I can think is stuff like identifying promoters or mutations in promoters, i.e. non-coding regions with a known function and can be trained against appropriate independent experimental observations, e.g. gene expression. If the k-mers were large it might work for RNA secondary structure because it would trap compensatory mutations. Vectorisation for individual nucleotides, i.e. dumbie variables, would not be successful at trapping compensatory mutations.
In summary it has niche application.
BERT is a famous large language model BTW.

If k-mer =3 that's about non-linearity of proteins, but it depends how the sliding window is deployed. It depends on the step size of the sliding window.
To try and answer this better, proteins fold back on themself triggering "co-evolution" or more accurately "covariance" between non-linear sites. This defines the triplet codon:

in machine learning this is needed because it will struggle to find it by itself
it will help it identify co-evolving mutations at non-linear locations across the gene

It might also identify selection events better.
Summary If you've sat down and thought the equation the RNN is attempting to unravel its mind-boggling. More worrisome does it really do it? It's really complicated to simply discuss, but when you  give it key information the ML/DL gets there quicker. The authors will be aware of that in their training and testing.
Honestly, we're at the limits of reasonable application of ML/DL without going into RasNet type infrastructure. My system is different the idea is the same - give the training model the key biological information, it gets to the answer quicker, a lot quicker.

I am sorry, but this does not answer my question regarding why one would use k-mers for k>1 versus 1-mers with positional encoding as is the norm with transformers. Not just BERT and other NLP transformer-like models (aka LLMs) but even ViT uses positional encoding

I've tried my best above.

[I'm not an expert in this field, this is only my "best try". Maybe you could get better mathematically motivated answers on ai.stackexchange]
First note that, as pointed in the DNA BERT paper, k-mers are a very common approach to working with sequences. This is true both for older approaches as well as newer Transformer-type models. For example the recent Nucleotide Transformer, DNAGPT or MuLan-Methyl also use k-mers.
While I'm not aware of a published systematic comparison between k-mers and single bases, I'm sure I've seen papers that vary k as a hyperparameter (so would notice if k=1 gives best results), and I have to believe someone somewhere at some point has thought about not using k-mers.
There are some Deep Learning papers that use individual one hot-encoded nucleotides, for example SpliceAI, but they typically only look at relatively short, fixed-length, sequences (relatively: this paper does go up to 10k).
As to the "why", my (limited) understanding is that there are two factors. If you think of NLP, a token is typically a word from the language. Taking each nucleotide as a token means you consider a language with only 4 words, that is a very limited amount of information. It does make sense to use more meaningful tokens: to respond to your comment, even ViT uses patches (and does not attempt to work on individual pixels). Based on our knowledge of biology, we can expect that most sequences have meaningful "motifs" that are made of a few consecutive nucleotides. If you feed the nucleotides to the Transformer, it has to use attention to detect these motifs as well as their long-range interaction. If instead you use a "patch" (or k-mer) that captures entire motifs, this can be learned by an initial linear projection and you only use the attention for long-range interactions.
Finally, the other reason is computing power, as hinted to in the DNAGPT paper:

In the non-overlapped k-mers strategy, the shift is equal to K, resulting in N/k tokens for an N-length sequence and improving the efficiency by k times.

Maybe taking 1-mers could work just as well as 6-mers, but the number of relationships to consider increases even more, a non-overlapping k-mer strategy does decrease the of interactions to consider, and an overlapping k-mer strategy where each patch is preprocessed (e.g. with a linear model) can also increase efficiency (so to speak the model does not need to consider short-range interactions which are already accounted for).

I would have actually given this an upvote if I hadn't just discovered DNABERT-2 where they demonstrate the exact opposite case; That k-mers are computationally inefficient. Also, I have a hard time understanding why this would be so different from proteomics. As I wrote to M_, "You don't need k-mers to use language model; If you tokenize 1-mers, you essentially have a language with four words and can thus use language models. I've done this in proteomics before (having a 21 word language, one for each amino acid, in that case)."

Since you bring proteins into the question: There are some differences between DNA and protein sequences with regards to their biology though. Protein sequences are in many cases the result of stringent selection, while for DNA sequences there is more sequence space that may not have "meaning". Large portions of the genome are non-coding, some of that sequence space is regulatory, but there is also a large fraction that is basically broken transposons. Was your attempt using a "1 amino acid - 1 word" approach in proteomics successful / better than other approaches?

@Niklas Are there other approaches? Even AlphaFold one-hot-encodes amino acids, which is the simplest "1 amino acid - 1 word"-approach out there (and more practical in proteomics than in NLP since there aren't as many amino acids as there are words in languages)

Sure there are other conceivable approaches, you could encode amino acids based on various chemical properties (charge, hydrophobicity, side-chain size). Especially for proteomics, you could also expand the alphabet to included chemical modifications. I do not know if they have been tried. ProteinBERT uses a 26-letter alphabet: https://academic.oup.com/bioinformatics/article/38/8/2102/6502274?login=false#394300229
But I feel we are getting a bit away from your initial question, which I think is difficult to know, unless you ask the creators of DNABert why they chose kmers.

For proteins, note also that even with one-hot encoding, you get an alphabet with 22 "words", which is already much better than 4. The DNABERT-2 paper also doesn't use one-hot encoding, but (see section 3.1 and Fig. 3) builds a vocabulary with 4096 "words" of length 6-7 nucleotides. So I think my tentative answer still stands: single nucleotide encoding (1-mers) hasn't appeared to be better in most of the literature. Note also (in their intro) their reported problem is that they used overlapping k-mers, which i)leaks information and ii) increases length

As for other tokenization approaches in DNA, you can also look at this, that, or that, with the caveat that I haven't seen these approaches used much in other papers.

Additional clarification: with the overlapping k-mers from DNABERT, you have as many k-mers as nucleotides in the sequence (minus k), so the "computational efficiency" argument of DNABERT-2 applies as well to 1-mers. What is more efficient is either non-overlapping k-mers or other tokenizations that create groups of nucleotides (as in DNABERT-2 or seq2vec).

@Alexlok Actually, I kind of regret how I have phrased my objections in the comments and the summary in the title; The body of the question is fine though. My objection does not apply to non-overlapping k-mers (since then we are not talking about using a rolling window on DNA sequences anymore; So, there's not redundancy in repeating information). I'll edit the title to make this clearer.


text:
Looping over a list of cell models in Python
I have two data sets

Mutations in 10 genes
Gene expression in 18 cancer 3D models

For every (mutation, gene) and gene expression pair we run an ANOVA test
This is my data

> mut=read.delim("mut.tsv",row.names = 1)
> expr=as.data.frame(t(read.delim("expr.tsv",row.names = 1)))

dput(expr)

structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723,
-0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448,
0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995,
-1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465,
-0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804
), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838,
-0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099,
-0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653,
-1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086,
-1.48513703581865, -1.48373348870346, -1.76484058374354), class = "data.frame", row.names = c("ModelM",
"ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH",
"ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU",
"ModelT2", "ModelC2", "ModelS2", "ModelX"))

dput(mut)

structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L,
0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L,
0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L,
0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L
), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L,
1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L,
1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L,
0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L,
0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L,
0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L,
0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L,
0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L,
0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),
    ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut",
"Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut",
"Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

To test this, for the first 3D cancer model as instance, I have done like below
           GeneA_KO   GeneB_KO
ModelM   0.07563831 -0.8701573
ModelO  -0.16284197 -0.3590087
ModelS  -0.11070310 -0.9904069
ModelA  -0.45715462 -0.1206570
ModelI  -0.50581519 -1.4142763
ModelC   0.36607278 -0.8664211
ModelT   0.28865046 -0.8169361
ModelH  -0.02876465 -0.8945972
ModelE  -0.69073282 -1.0008458
ModelR   0.15736612 -0.8507545
ModelA2 -1.29451547 -1.0531424
ModelP   0.12416743 -1.0758068
ModelE2 -0.16804531 -1.2470832
ModelU   0.01865049 -1.0409846
ModelT2 -0.49588083 -0.2323887
ModelC2 -0.20666863 -1.4851370
ModelS2 -0.14860949 -1.4837335
ModelX   0.06684453 -1.7648406

           ModelM ModelO ModelS ModelA ModelI ModelC ModelT ModelH ModelE
Gene1_mut       0      0      0      0      0      0      0      1      1
Gene2_mut       0      0      0      0      0      0      0      1      0
Gene3_mut       0      0      0      0      0      0      0      1      0
Gene4_mut       1      1      1      1      1      1      1      1      0
Gene5_mut       0      0      0      0      0      0      0      0      1
Gene6_mut       0      0      0      0      0      0      1      1      0
Gene7_mut       0      0      0      0      0      0      0      1      0
Gene8_mut       0      0      0      0      0      0      1      1      0
Gene9_mut       0      0      0      0      0      0      1      0      0
Gene10_mut      0      0      0      0      0      0      0      0      1
           ModelR ModelA2 ModelP ModelE2 ModelU ModelT2 ModelI2 ModelC2 ModelS2
Gene1_mut       1       0      0       1      1       0       1       1       0
Gene2_mut       0       0      0       0      0       0       0       1       1
Gene3_mut       0       0      0       0      1       0       0       0       0
Gene4_mut       1       1      1       1      1       1       1       0       1
Gene5_mut       0       0      0       0      0       0       0       0       1
Gene6_mut       0       0      0       1      1       0       0       1       1
Gene7_mut       1       0      0       0      0       0       0       0       1
Gene8_mut       0       0      0       0      1       0       0       1       0
Gene9_mut       1       0      0       0      1       0       0       0       0
Gene10_mut      0       0      1       0      1       0       0       0       1

> my_data=cbind(mut[,1],expr[,1])
> colnames(my_data)=c("mut","expr")
> res.aov <- aov(expr ~ mut, data = as.data.frame(my_data))
> # Summary of the analysis
> summary(res.aov)
            Df Sum Sq Mean Sq F value Pr(>F)
mut          1  0.145  0.1448   0.192  0.673
Residuals    8  6.022  0.7527

Now I do a for loop for all of the 18 cancer cell lines
> for(i in 1:ncol(g))
+ {
+   my_data=as.data.frame(cbind(m[,1],g[,1]))
+   #tidy will summarise and return neat format
+   print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
+ }
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA
# A tibble: 2 × 6
  term            df sumsq meansq statistic p.value
  <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
1 my_data[, 1]     1 0.145  0.145     0.192   0.673
2 Residuals        8 6.02   0.753    NA      NA

I done this in R but of course I got when doing in Python
Python
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
os.getcwd()
os.chdir("/Users//Downloads/test")
os.getcwd()
m = pd.read_csv("mut.tsv", sep="\t", index_col=0,delim_whitespace=True)
m=pd.DataFrame(m)
m
g = pd.read_csv("expr.tsv", sep="\t",index_col=0,delim_whitespace=True).T
g=pd.DataFrame(g)
g
lst = list(g.columns)

lst
dfs = []
for Model in lst:
   dfs.append(pd.concat([m.iloc[:, 0].reset_index(drop=True), g.iloc[0,:].reset_index(drop=True)],axis=1))
   dfs.append(df.dropna())
   dfs.append(df.columns==['mut','expr'])
   model = ols('expr ~ C(mut)', data=df).fit()
   aov_table = sm.stats.anova_lm(model, typ=2)
   aov_table
   final = pd.concat(dfs, ignore_index=True)

    for Model in lst:
      Input In [165]
        for Model in lst:
                         ^
    IndentationError: expected an indented block

When asking about errors, please remember to always include the actual error. We cannot help you debug errors you don't share.

thank you, I have added the error in my main post

@Angel, you still did not include the error that you get when you are using R. Anyway, why would you define i in your loop and not use it? And when it comes to the Python error, it is obvious!

Thank you @haci I get error in Python part only. The R code is OK

Thanks I have edited my post so I believe the R code will run smoothly. I have added the desire output from R code to the post. First, I have tried for one sample then using a for loop to doing one-way ANOVA over 18 samples. But, I do not know how to do this by Python

Here is one potential approach:
Example data:
expr1 <- structure(list(GeneA_KO = c(0.0756383119639074, -0.162841965682723, -0.110703100204517, -0.457154616239399, -0.505815193885352, 0.366072775033448, 0.288650458935083, -0.0287646498640309, -0.690732824311905, 0.157366122620995, -1.29451546737684, 0.124167425460344, -0.168045313861751, 0.0186504943306465, -0.495880827400706, -0.206668627339351, -0.148609490335447, 0.066844529089804), GeneB_KO = c(-0.870157285953621, -0.35900872537626, -0.990406888128838, -0.120656981579292, -1.41427631791049, -0.866421095836781, -0.816936112146099, -0.894597236248507, -1.00084584207823, -0.850754488133456, -1.05314239761653, -1.07580675428372, -1.24708324665688, -1.04098460460899, -0.232388675985086, -1.48513703581865, -1.48373348870346, -1.76484058374354)), class = "data.frame", row.names = c("ModelM", "ModelO", "ModelS", "ModelA", "ModelI", "ModelC", "ModelT", "ModelH", "ModelE", "ModelR", "ModelA2", "ModelP", "ModelE2", "ModelU", "ModelT2", "ModelC2", "ModelS2", "ModelX"))
mut <- structure(list(ModelM = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelO = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelS = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelA = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelT = c(0L, 0L, 0L, 1L, 0L, 1L, 0L, 1L, 1L, 0L), ModelH = c(1L, 1L, 1L, 1L, 0L, 1L, 1L, 1L, 0L, 0L), ModelE = c(1L, 0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 1L), ModelR = c(1L, 0L, 0L, 1L, 0L, 0L, 1L, 0L, 1L, 0L), ModelA2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelP = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 1L), ModelE2 = c(1L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), ModelU = c(1L, 0L, 1L, 1L, 0L, 1L, 0L, 1L, 1L, 1L), ModelT2 = c(0L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelI2 = c(1L, 0L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L), ModelC2 = c(1L, 1L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L),     ModelS2 = c(0L, 1L, 0L, 1L, 1L, 1L, 1L, 0L, 0L, 1L)), class = "data.frame", row.names = c("Gene1_mut", "Gene2_mut", "Gene3_mut", "Gene4_mut", "Gene5_mut", "Gene6_mut", "Gene7_mut", "Gene8_mut", "Gene9_mut", "Gene10_mut"))

R code:
my_data=cbind(mut[,1],expr1[1:10,1])
colnames(my_data)=c("mut","expr1")
res.aov <- aov(expr1 ~ mut, data = as.data.frame(my_data))
summary(res.aov)
#>             Df Sum Sq Mean Sq F value Pr(>F)
#> mut          1 0.1364  0.1364   1.114  0.322
#> Residuals    8 0.9789  0.1224

for (j in 1:nrow(mut)) {
  for(i in 1:ncol(expr1)) {
    my_data=as.data.frame(cbind(mut[,j],expr1[1:10,i]))
    print(broom::tidy(aov(my_data[,2] ~ my_data[,1], data = my_data)))
  }
}
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.136  0.136      1.11   0.322
#> 2 Residuals        8 0.979  0.122     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.541 0.541       7.38  0.0264
#> 2 Residuals        8 0.587 0.0733     NA    NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0612 0.0612     0.465   0.515
#> 2 Residuals        8 1.05   0.132     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0637 0.0637     0.479   0.508
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.246  0.246      2.26   0.171
#> 2 Residuals        8 0.869  0.109     NA     NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.313  0.313      3.07   0.118
#> 2 Residuals        8 0.815  0.102     NA     NA
#> # A tibble: 2 × 6
#>   term            df   sumsq  meansq statistic p.value
#>   <chr>        <dbl>   <dbl>   <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.00108 0.00108   0.00777   0.932
#> 2 Residuals        8 1.11    0.139    NA        NA
#> # A tibble: 2 × 6
#>   term            df sumsq meansq statistic p.value
#>   <chr>        <dbl> <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.220  0.220      1.94   0.201
#> 2 Residuals        8 0.907  0.113     NA     NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0529 0.0529     0.398   0.546
#> 2 Residuals        8 1.06   0.133     NA      NA
#> # A tibble: 2 × 6
#>   term            df  sumsq meansq statistic p.value
#>   <chr>        <dbl>  <dbl>  <dbl>     <dbl>   <dbl>
#> 1 my_data[, 1]     1 0.0901 0.0901     0.695   0.429
#> 2 Residuals        8 1.04   0.130     NA      NA

# write.csv(mut, "mut.csv", quote = FALSE, row.names = FALSE)
# write.csv(expr1[1:10,], "expr1.csv", quote = FALSE, row.names = FALSE)

Created on 2023-09-19 with reprex v2.0.2

Python code:
import os
import pandas as pd
from statsmodels.formula.api import ols
import statsmodels.api as sm
m = pd.read_csv("mut.csv", sep=",")
m=pd.DataFrame(m)
m.info()
g = pd.read_csv("expr1.csv", sep=",")
g=pd.DataFrame(g)

lst = list(g.columns)
m_lst = list(m.columns)

dfs = pd.DataFrame({})
for Mut in m_lst:
    for Model in lst:
        df = pd.concat([m[[Mut]].reset_index(drop=True), g[[Model]].reset_index(drop=True)],axis=1)
        df = df.rename(columns={Mut: "mut", Model: "g"})
        model = ols('g ~ C(mut)', data=df).fit()
        aov_table = sm.stats.anova_lm(model, typ=2)
        aov_table[['Model']] = Model
        aov_table[['Mut']] = Mut
        dfs = dfs.append(aov_table)

with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    print(dfs)
            sum_sq   df         F    PR(>F)     Model      Mut
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelM
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelM
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelM
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelM
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelO
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelO
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelO
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelO
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelS
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelS
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelS
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelS
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelA
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelA
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelA
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelA
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelI
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelI
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelI
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelI
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO   ModelC
Residual  0.978877  8.0       NaN       NaN  GeneA_KO   ModelC
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO   ModelC
Residual  0.586718  8.0       NaN       NaN  GeneB_KO   ModelC
C(mut)    0.061205  1.0  0.464539  0.514753  GeneA_KO   ModelT
Residual  1.054037  8.0       NaN       NaN  GeneA_KO   ModelT
C(mut)    0.063734  1.0  0.479234  0.508361  GeneB_KO   ModelT
Residual  1.063933  8.0       NaN       NaN  GeneB_KO   ModelT
C(mut)    0.245964  1.0  2.263619  0.170857  GeneA_KO   ModelH
Residual  0.869278  8.0       NaN       NaN  GeneA_KO   ModelH
C(mut)    0.312937  1.0  3.072788  0.117706  GeneB_KO   ModelH
Residual  0.814730  8.0       NaN       NaN  GeneB_KO   ModelH
C(mut)    0.001082  1.0  0.007771  0.931920  GeneA_KO   ModelE
Residual  1.114159  8.0       NaN       NaN  GeneA_KO   ModelE
C(mut)    0.220171  1.0  1.940909  0.201060  GeneB_KO   ModelE
Residual  0.907496  8.0       NaN       NaN  GeneB_KO   ModelE
C(mut)    0.052891  1.0  0.398296  0.545575  GeneA_KO   ModelR
Residual  1.062350  8.0       NaN       NaN  GeneA_KO   ModelR
C(mut)    0.090105  1.0  0.694741  0.428739  GeneB_KO   ModelR
Residual  1.037562  8.0       NaN       NaN  GeneB_KO   ModelR
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelA2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelA2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelA2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelA2
C(mut)    0.004637  1.0  0.033399  0.859537  GeneA_KO   ModelP
Residual  1.110605  8.0       NaN       NaN  GeneA_KO   ModelP
C(mut)    0.276724  1.0  2.601573  0.145422  GeneB_KO   ModelP
Residual  0.850943  8.0       NaN       NaN  GeneB_KO   ModelP
C(mut)    0.044310  1.0  0.331000  0.580888  GeneA_KO  ModelE2
Residual  1.070932  8.0       NaN       NaN  GeneA_KO  ModelE2
C(mut)    0.170278  1.0  1.422852  0.267110  GeneB_KO  ModelE2
Residual  0.957389  8.0       NaN       NaN  GeneB_KO  ModelE2
C(mut)    0.001687  1.0  0.012120  0.915049  GeneA_KO   ModelU
Residual  1.113555  8.0       NaN       NaN  GeneA_KO   ModelU
C(mut)    0.008679  1.0  0.062048  0.809563  GeneB_KO   ModelU
Residual  1.118988  8.0       NaN       NaN  GeneB_KO   ModelU
C(mut)    0.136365  1.0  1.114460  0.321950  GeneA_KO  ModelT2
Residual  0.978877  8.0       NaN       NaN  GeneA_KO  ModelT2
C(mut)    0.540949  1.0  7.375925  0.026417  GeneB_KO  ModelT2
Residual  0.586718  8.0       NaN       NaN  GeneB_KO  ModelT2
C(mut)    0.017610  1.0  0.128353  0.729416  GeneA_KO  ModelI2
Residual  1.097631  8.0       NaN       NaN  GeneA_KO  ModelI2
C(mut)    0.260821  1.0  2.407079  0.159387  GeneB_KO  ModelI2
Residual  0.866846  8.0       NaN       NaN  GeneB_KO  ModelI2
C(mut)    0.191207  1.0  1.655405  0.234206  GeneA_KO  ModelC2
Residual  0.924035  8.0       NaN       NaN  GeneA_KO  ModelC2
C(mut)    0.033474  1.0  0.244741  0.634101  GeneB_KO  ModelC2
Residual  1.094193  8.0       NaN       NaN  GeneB_KO  ModelC2
C(mut)    0.044621  1.0  0.333425  0.579533  GeneA_KO  ModelS2
Residual  1.070620  8.0       NaN       NaN  GeneA_KO  ModelS2
C(mut)    0.096956  1.0  0.752533  0.410940  GeneB_KO  ModelS2
Residual  1.030711  8.0       NaN       NaN  GeneB_KO  ModelS2
```

Sorry I was running your python code but I get this error for Mut in m_lst:
Input In [81]
for Mut in m_lst:
^
IndentationError: expected an indented block

The IndentationError means you haven't got the correct whitespace - you need 4 spaces in front of each indented line @Angel and no whitespace in front of 'non-indented' lines (basically, copy it exactly as I have it in my answer and it will work)

Thanks a lot that's weird I just copied and pasted that but gives error

You're using spyder, so perhaps try going to "Source --> Fix Indentation" and see if your code changes and if the error persists

Source fixed Indentation error now I get this Traceback (most recent call last):
File "/Users/angel/Downloads/test/test.py", line 20, in
model = ols('g ~ C(mut)', data=df).fit()
File "/Users/angel/Downloads/docker/.venv/lib/python3.9/site-packages/statsmodels/base/model.py", line 209, in from_formula
raise ValueError('endog has evaluated to an array with multiple '
ValueError: endog has evaluated to an array with multiple columns that has shape (10, 18). This occurs when the variable converted to endog is non-numeric (e.g., bool or str).

That error is harder to understand, sorry @Angel; I think you should post another question with a reproducible example showing how you got that error so we can figure out what the problem is.

Before using this-->  with pd.option_context('display.max_rows', None, 'display.max_columns', None):,  print(dfs) could display the 1st 10 rows,  tried, worked!

Moreover, after using this--> python -V,  my z shell on my macbook air is using Python 3.9.16 .

Apparently the pd.option_context('display.max_rows', None, 'display.max_columns', None) approach doesn't work with jupyter notebooks; see https://stackoverflow.com/questions/16424493/pandas-setting-no-of-max-rows for more details @envs_h_gang_5


text:
Metrics to use to estimate variant library evenness, uniformity, or bias
I have a 15k library of variants (OTUs) from a genome region. I want to compare the evenness of this library and understand of the library bias.
I tried different methods to calculate evenness or uniformity of this library such as Shannon evenness index, relative abundance rank plot, kurtosis , skewness. I am confused with the results and why they are not aligned. For example my Shannon evenness index is 0.98 which means library is very even, while my kurtosis is very large (120), indication not uniform distribution.
I should note that for calculating kurtosis I use relative abundance which is read count of every variant /total read count of all variants in that library. My question is what metrics I can use to have a better understanding of library evenness or uniformity, other that above metrics, and why my results are so contradictory?

Shannon and Kurtosis are very different statistics. Kurtosis is generally used in conjunction with skewness as a very general approach to normal distribution signal. I would simply say they are measuring different things


text:
Export multiple variables from shell script into a Nextflow output channel
Within a Nextflow workflow, some process runs a shell script that creates multiple variables, e.g.
var_a=$(echo "1")
var_b=$(echo "2")
var_c=$(echo "3")
var_d=$(echo "4")

These variables contain statistics (e.g. number of input reads, quality scores, ...) that I want to keep until the end of the pipeline, to then save them and perform QC. So I'd like to export them out of this process and carry them to the end of the workflow. To be clear, the number of variables is not unknown, it's just too many for convenience if I need to list them all every time.
I understand that I can use the env qualifier or an intermediate file to export the variables individually:
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d)

  shell:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow{
  get_vars | view
}

with output:
[ea/a145b4] process > get_vars [100%] 1 of 1 ✔
[1, 2, 3, 4]

However that becomes very impractical if I have a lot of such variables (that I do want to continue passing along to later processes, along with a much longer tuple of other relevant variables).
My question: is there a convenient way to "encapsulate" all these variables in the output channel, for example with a Groovy map similar to the meta map in nf-core? So I could directly write something like:
process get_vars {
  output:
    val(my_vars)

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''

  my_vars = [:]
  my_vars.a = var_a
  my_vars.b = var_b
  my_vars.c = var_c
  my_vars.d = var_d
}

workflow{
  get_vars | view
}

which obviously fails because var_a is not defined outside of the ''' block.

Can't you simple write all vars to a plain text file and then have some groovy/nf code to loop over the content and pass the it to any data structure that is appropriate, be it a map, list, anything?

That does sound like a solution, at least I would only have to pass along a path. But that means parsing again a file of variables I have already parsed, which doesn't feel like the most direct approach. If you post it as an answer I would accept it.

I think you can use output for this, ouput: env varname, as is mentioned here https://www.nextflow.io/docs/latest/process.html#outputs

I am not sure if i get what you mean by encapsulate, you can give your tuple a name using emit?
process get_vars {
  output:
    tuple env(var_a), env(var_b), env(var_c), env(var_d), emit: vars

  script:
  '''
  var_a=$(echo "1")
  var_b=$(echo "2")
  var_c=$(echo "3")
  var_d=$(echo "4")
  '''
}

workflow {
  get_vars()
  get_vars.out.vars.view()
}

That gives:
 [1, 2, 3, 4]

If I'm not mistaken, using emit (described in the DSL2 docs) is to create a separate (named) channel, so I still have to ensure that it keeps the order in sync with other channels, right?

Here is an illustration keeping in sync with a main channel:
`process get_vars {
input:
val(main_chan)
output:
tuple env(var_a), env(var_b), env(var_c), env(var_d), emit: vars
val(main_chan), emit: main_chan
shell:
'''
var_a=!{main_chan}
var_b=$((!{main_chan} + 1))
var_c=$((!{main_chan} + 2))
var_d=$((!{main_chan} + 3))
'''
}`
`process next_process {
input:
val(vars)
val(main_chan)
output:
tuple(val(vars), val(main_chan))
script:
'''
'''
}`

workflow {   main_chan = Channel.of( 1, 2, 3 )   get_vars(main_chan)   next_process(get_vars.out.vars, get_vars.out.main_chan)   next_process.out.view() }

The key here is that next_process(get_vars.out.vars, get_vars.out.main_chan) needs to receive both channels from the previous process at once. Feel free to include this in your answer if you think it would be relevant

I see, you can use a id (like sample-name or meta.id in nfcore) that you include in both channels and then use to join them together into a single channel to make sure that information for one sample moves though the pipeline together.


text:
Changing transparency of markers in Chimera X
I am using the software Chimera X to visualize the 3d composition of molecules. I am very new to Chimera and only use it for visualization. The files I'm using are cmm files (xml-formatted) basically looking like below:
...
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
...
<link id1="0" id2="1" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="1" id2="2" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="2" id2="3" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="3" id2="4" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="4" id2="5" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
...

i.e. just a simple polymer consisting of alot of markers (spheres) and links between them (cylinders), and with colors for each one specified in rgb (see here). This is all well and good, but now I would like to set transparency for a subset of the markers and associated links. Preferably I would set this directly in the attributes in the cmm files, but I cannot find how to do this. The manual does mention that colours in general can be set with rgba where a means alpha level, but when trying to set a like seen in the xml block above, nothing happens (I tried values close to both 0 and 1; no difference)
The manual documents a function called 'transparency' which I tried to use with for example transparency 50,a #0:10-1200
to make atoms 10-1200 50% transparent. But I just get the error: "Missing or invalid "percent" argument: Expected a number". If I only type transparency 50 everything does turn transparent, however, I wouldn't really call it transparent since I cannot really see through them to see the other molecules on the other side, they just turn bleak.
Questions:

(Would someone explain how I could simply change the actual transparency of a large set of markers and links) (preferably by loading some file with specified marker id:s)? Also, is it possible to assign this directly in the cmm files, and if so, how?

UPDATE: The command transparency #1:2492-end 90 'All' worked to set the transparency of markers e.g. 2492-end. What I would still like to know is if it is possible to set this directly in the cmm-files.

RGBA color values are an extension of RGB color values with an alpha channel - which specifies the opacity for a color. An RGBA color value is specified with: rgba(red, green, blue, alpha). The alpha parameter is a number between 0.0 (fully transparent) and 1.0 (fully opaque). , what does the  a="0"  do in your xml formatted markers and links ?

So, the 'a="0"' was a naive attempt by me to set the alpha-level directly in the xml file, which as I noted above did not work. I have since figured out how to set the transparency and will update my question.

you should try to report the bug a = 0 not working to chimerax developers. To report a problem with ChimeraX:
If possible, please use Help... Report a Bug in the ChimeraX menu. Otherwise (for example, if ChimeraX will not start), send e-mail to chimerax-bugs@cgl.ucsf.edu, making sure to include any information needed to reproduce the problem: program version, specific data, actions/steps taken, etc. The information will be publicly accessible, so do not include data or descriptions that you wish to remain private.

from https://www.cgl.ucsf.edu/chimerax/docs/contact.html

I will look in to that, however I am not certain it is a bug: Even though it seems very reasonable to be able to set it that way, that was a guess on my part. The general support of rgba colors is however otherwise documented in the manual. Regarding your suggested edit: The markers are literally referred to as 'atoms' in Chimera, and it is those which I wish to make transparent. (I am however changing it to 'marker' as a compromise =).)

A marker is a pseudoatom, typically used to mark the position of some feature within volume data; markers are implemented as atoms, so actions that apply to atoms in ChimeraX also apply to markers: they can be specified in commands, selected, etc. https://www.cgl.ucsf.edu/chimerax/docs/user/markers.html

dont know why I didnt like ChimeraX at the time , ages ago,  but haved worked only with just plain Chimera, had to dig an old VMwareVM to take these few pics

users mailing list is https://mail.cgl.ucsf.edu/mailman/mailman3/lists/chimera-users.cgl.ucsf.edu/ .. dont get why  a='50' or a='0.5' dont work with chimera too

Ok, this is not an answer, but to better visualize the problem I recreated the question/answer on Chimera that has a different syntax for the transparency  command see : https://www.cgl.ucsf.edu/chimera/docs/UsersGuide/midas/transparency.html and https://www.cgl.ucsf.edu/chimerax/docs/user/commands/transparency.html .
Input file is marker.cmm , to be opened in FILE MENU of VOLUME TRACER from VOLUME MENU :
<marker_set name="marker set 1">
<marker id="14" x="0.018" y="-0.180" z="0.215" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="15" x="-0.032" y="-0.191" z="0.241" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="16" x="-0.033" y="-0.117" z="0.246" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="17" x="0.020" y="-0.138" z="0.244" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="18" x="0.047" y="-0.083" z="0.280" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="19" x="0.073" y="-0.061" z="0.240" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<marker id="20" x="0.010" y="-0.068" z="0.271" r="1.000000" g="0.254900" b="0.050980" a="0" radius="0.018754" />
<link id1="14" id2="15" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="15" id2="16" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="17" id2="18" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="18" id2="19" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
<link id1="19" id2="20" r="0.800000" g="0.800000" b="0.000000" a="0" radius="0.006251"/>
</marker_set>

Image before command :

Image after transparency  50,a,b #0:16-end
(#0 being 0 model as the only loaded model on my Chimera instance) :

Better transparency visualizations :

Dont get why a='50' or a='0.5' in marker,link definition dont work with chimera too

Together with my edited question I think this all suffice as "good enough" for now at least. I am satisfied for the moment just using something like transparency #1:2492-end 90 'All', but I guess when I get more complicated regions it would be good to be able to load the marker spec, i.e. e.g. #1:2492-2500, 3000-3500, ... from a file if that is possible. But I guess one can run command-files right? and I suppose there's no limit on the length on the spec, so perhaps I could just create one very long command in a file. (Many thanks for your efforts!!! unfortunately I need more rep to upvote)