text:
How to visualize unique and overlapping genes in upSet R plot
I am still a newbie in the field of bioinformatics. I have been attempting to visualize unique and common genes across six different tissues using the UpSetR plot package. This script was used to generate the plot; however, out of the six tissues, only six could be on the plot.
library(data.table)
library(UpSetR)
d <- data.table::fread(text = "Brainstem    Cerebellum  Olfactory Bulb  Heart   Kidney  Liver
Rac3    Coq9    Ndufs7  Slc25a10    Aacs    Acss3
Ndufb5  Mtnd5   Cotl1   Adhfe1  Aass    Cyp4v2
Clic6   Kpna4   Gabarapl2   Ddx5    Abcc4   Slc7a2
Rbmx    Snrpg   Rbm3    Pgd Acot12  Ehhadh
Ptcd3   Micos13 Atp5f1e Gss Acsf2   Selenbp2
Ndufa13 Ndufa5  Ndufb3  Arl3    Acsm3   Ugt2b17
Amt Ndufa10 Cstf3   Shmt2   Acy3    Gstp1
Ndufs1  Bcr Mtco2   Tpt1    Adk Pdk1
Acan    Mpc1    Bcap31  Me1 Ahcyl1  Bgn
Wasl    Cfdp1   Ndufaf2 Abcg2   Aif1l   Col4a2
Rbm3    Rps15   Gfap    H1-1    Akr1b1  Ttn
Mtnd5   Snu13   Atp6v0c Pcbp2   Amacr   Cox6b1
Ddc Lin7c   Tmed10  Ddx17   Atp2b1  Mup2
Tmc7    Ndufb5  S100a5  Slc25a1 Ca5b    Cyp2d9
Erc1    Ampd3   Ap3b1   H1-3    Cbr1    Cyp2b10
Nup35   Ipo5    Gja1    Cstb    Cbr3    Ftl1
Ndufb2  Ndufb10 Ndufb7  Tst Ces1f   Cyp8b1
Ndufb1  Ebp Ciapin1 Vat1    Cnpy2   Rbm3
Rexo2   Rabac1  Ndufb4  Glul    Col1a1  Ehd1
Cbln3   Wdfy2   Adk Rpl29   Col1a2  Cyp4a10
Arfgef3 Clcn6   Gfpt1   Pc  Col4a1  Ndufaf2
Dagla   Vbp1    Rpl28   Acy3    Col4a2  Cops2
Slc25a31    Nae1    Dgkz    Akr1a1  Cox7a1  Cyp2a4
Tph2    Gucy1a1 Tpbg    Pecr    Cpne5   Rbp1
Actbl2  Ap1s1   Gprc5b  Hnrnpd  Creg1   S100a13
Ndufb9  Ndufs1  Acyp2   Serpina1c   Crot    Gsta4
Timm10b Ndufa3  Kirrel2 Hsd17b11    Cryab   Ldhb
Syt12   Ndufs4  Ndufb11 Rpn1    Cyp2d26 Dhtkd1
Unc13c  Dnajb4  Ndufb10 Acly    Cyp2d9  Ndufb7
Serpina1d   Chmp4b  Wdr6    Atp6v1b2    Cyp2j6  Acot13
Yme1l1  Pdcd6   Mpst    Pxmp2   Cyp4a10 Fmo3
Aprt    Ndufaf2 Smdt1   Arpc3   Cyp4a14 Clta
Mtnd1   Igsf21  Sec61g  Kpna6   Cyp4b1  Eftud2
Jpt2    Stxbp3  C1qb    Atp6v1g1    Cyp7b1  Tex15
Ndufv3  Ndufb1  Ndufa2  Slc25a31    Ddc Ndufs1
Chordc1 Tpd52l2 Hddc2   Aldh3a2 Dpyd    Sult1c2
Cbln1   Sh3gl1  Agfg1   Gng12   Echdc3  Srl
Rbbp4   Pnn Manf    Tmpo    Eea1    Aimp1
Ndufaf2 Prodh   Kcnj10  Acox1   Ehhadh  Cyp2b9
Egfr    Exoc7   Ptgds   Rhoa    Folh1   Cyp4a14
Ca4 Ubxn1   Rpa3    Mapre1  Fth1    Csad
Vps33b  Nfu1    Eif3e   Gk  Ftl1    Gck
Ndufb11 Etnppl  Isg15   Rps2    Gsta2   Ss18l2
Gatad2b Spock1  Tmsb10  Atp6v1e1    Gsta3   Hebp1
Scn8a   Tagln2  Chmp1a  Sult1d1 Gsta4   Mup20
Ensa    Cbln1   Mobp    Plin1   Hacl1   Nudt7
Hkdc1   Sart3   Mtmr9   Idh1    Hmgb2   Tmem19
Ndufa5  Ndufs6  Ptprf   Cpne1   Hmgcs2  Mup18
Prpf6   Tagln   Hmox2   Prodh   Hmgn2   Selenbp1
Glipr2  Ndufs8  Ptprd   Asl Hpd Mpp1
Serpina3k   Drc7    Pmm1    Stmn1   Hsd3b2  Gstt2
Serpina1a   Hspg2   Pde6d   Abhd14b Iqgap2  Tmco1
Camk2g  Aamdc   Ipo7    Cryz    Krt18   Slc30a5
Nmnat1  Rpl26   Ndufc2  Amacr   Krt7    Hsd17b6
Tbcd    Nmral1  Csde1   Hmgb2   Krt8    Dusp3
Ckm Uck1    Itpa    Sfxn1   Ldhd    Aass
Mtif2   Lysmd2  Ndufa5  Acad11  Lmnb2   Sh3bgrl3
Hmgn5   Sars2   Actr3b  Alox12  Lsm6    Tex44
Tmem245 Cep170b Meis2   H2ac20  Mapre1  Cyp7b1
Cbx3    Bzw1    Bcas3   Csad    Mbp Tymp
Gjc3    Apex2   Ensa    Ucp1    Mpv17l  Fth1
Ugt8    Nucks1  Serpina1a   Gpd1    Mrpl22  Serpina3m
Lsmem1  Ndufb6  Serpina3k   Xab2    Mrpl23  Ces3a
Rbsn    Cp  Tnks1bp1    Adh1    Mrpl40  Trmt112
Cpped1  Npm1    Prnp    Dcps    Mrps26  Rdh16
Arpp19  Rhoa    Uap1l1  Dhx9    Myh6    Ndufs4
Ndufa6  Ndufa8  Ndufs3  Ezr Ndufa10 Hspa1b
Ndufs2  Ndufb2  Ndufv2  Sri Ndufa12 Ndufv1
Mug1    Ndufv2  Ndufs4  Ptrh2   Ndufa2  Sult1d1
Ndufs8  Nek9    Cltb    Mecr    Ndufa5  Acot3
Ndufb7  Ndufa12 Ndufa12 Sqor    Ndufa6  Ndufa12
Ndufv1  Ndufs5  Wdfy2   Echdc2  Ndufa7  Upb1
Ndufs3  Ndufa2  Ndufa10 Fabp1   Ndufa9  Ttc39c
Dnah8   Lin7b   Mtnd1   Atp2a3  Ndufaf2 Sult1a1
Ndufc2  Nek7    Ctsa    Gsr Ndufb7  Sqor
Ndufa10 Ndufs3  Ndufv3  Adrm1   Ndufs1  Aox3
Ndufs7  Ndufv1  Ndufs6  Egfr    Ndufs2  Sult2a8
Ndufs5  Ndufb4  Clybl   Gmpr    Ndufs3  Hyou1
Gcfc2   Ndufa9  Stat1   Gdap2   Ndufs4  Egfr
Ndufb3  Ndufa13 Rpl23   Ywhaq   Ndufs6  Ndufa7
Ndufb4  Ndufv3  Ndufv1  Serpina1e   Ndufv1  Ndufa6
Ndufa12 Ndufs7  Rif1    Lta4h   Ndufv2  Ndufa9
Ndufb8  Ndufb3  Ndufb8  Kifbp   Ndufv3  Srd5a1
Serpina1e   Pgam2   Capg    Psma4   Nqo1    Mpc1
Ndufb10 Ndufc2  Amt Mtr Nucks1  Keg1
Ndufs6  Ndufa6  Ndufs1  Cps1    Nudt19  Serpina3k
Fam136a Cltb    Cep170b Ehhadh  Pals2   Gsta2
Ndufv2  Ndufb11 Ndufa13 Flii    Pecr    Gstt3
Ndufa11 Pcdhga4 Vwa1    Aldob   Ppid    Metap1d
Ndufa1  Mif Tinagl1 H1-4    Ptprk   Igf2bp3
Ndufb6  Ndufb9  Ndufa11 Kyat3   Pyroxd2 Ndufs5
S100b   H1-0    Actl6b  Coasy   Rbm3    Nudt19
Acta1   Rpl29   Ndufa8  Snrpe   Reep6   Atp6v1g1
Ndufa2  Ndufs2  Lsmem1  Vta1    Retsat  Maoa
Nfia    Acad8   Ilk Osbpl8  Sephs2  Mgst3
Ndufa7  Ndufa11 Ndufa7  Agpat2  Serpina3m   Comt
Ndufa8  Ech1    Ndufs5  Clic1   Sh3bgrl3    Mat1a
Ndufs4  Chmp2a  Msn Macroh2a1   Slc12a3 Acot4
Alas1   Micos10 Nup85   Sars1   Slc22a27    Fuom
Phyhd1  S100b   Nucks1  Pabpc1  Slc22a6 Mrpl14
Mtnd4   Synpr   Tceal3  H2ac21  Slc6a19 Serpinb1a
Ndufa9  Ndufa7  Ndufa9  Nudt8   Stat1   Ces3b
Cdh10   Gnpda2  Cep152  Krt18   Tor1aip1    Prelp
Ring1   Rbfox2  Etnppl  Fuom    Ugt1a2  Mup3
Tmsb10  S1pr1   Acta1   Cyb5b       Cyp2c37
Mrpl43  mt-Nd2  Myh11   Slc25a35        Col12a1
Icmt    Rnps1   Ddx46   Rbm3        Endov
Ncln    Ap3s1   Slc44a2 Gpd1l       S100a9
Gh1 Rbx1    Ndufa6  Lrrc59      Slc25a4
Bola3   Cat Ifit1   Hmgn2       Lum
Ahsg    Atp5if1 Vps52   Tkt     Myl1
Pdlim1  Rpl18a  Rac2    Hmgb1       Col1a1
Ca3 Cyp51a1 Tmco5b  H1-2        Tpm2
Crim1   Pkn1    Bclaf1  H4c1        Myh7
Supt4h1a    Ppp3r1  Ube2v2  Stat1       Dcn
S100a9  Sltm    Atp5if1 Aldh9a1     Cyfip1
Shisa6  Chmp7   Mcm6    Elavl1      Phf11
Adcy5   Bcl7c   Vac14   H2bc7       Cse1l
Slc4a2  Lama5   Aqp1    H1-5        Col1a2
Rgs8    Prnp    Hnrnpul1    Fah     Fabp3
Asph    Nup85   Dnajc5  Aldh1l2     Acta1
Hsd17b2 Iah1    Fam136a H3-5        Myh8
Prx Mbnl2   Ublcp1  H2ac4       Ogn
Bzw2    Rbm3    Cmtm5   Ldhd        Soat2
Mboat4  Hdgfl3  Nipa1   Cryl1       Sf3b5
Cep76   Hacd2   Bzw1    Ncl     Btk
Cacna1e Asph    Tgm2    Gpd2        Cela1
Ldhd    Lpcat4  Akr1e2  Ndrg1       Uap1
    Iqsec3  Nek7    Acacb       Actn2
    Kyat3   Gabbr1  Rcc2        Tnnt3
    Hnrnpul1    Pld2    Fgb     Upp2
    Atp5f1e C1qa    Ces1f       Tgtp2
    Mrps21  Irgm1   Gclc        Tpm1
    Itpa    Try5    Pvalb       Myh1
    Sec61g  Alg2    Bdh1        Tex264
    Acads   Pik3c3  Sord        Amy2a5
    Poldip2 Dynll2  Ssb     Cela2a
    Rpl36a  Ckm Fabp4       Ckmt2
    Slc6a5  Gstm5   Apoc4       Fmod
    Fn1 Urod    Bpnt1       Des
    Usp19   Rbbp7   Fasn        Ckm
    Prtn3   Ociad1  Hcfc1       Myl3
    Sh3bgrl Rpl29   Cyp4b1      Try5
    Gclm    Ngef    Ganab       Nrde2
    Egfr    Cryab   Tkfc        Ctrb1
    Mobp    Rnps1   Acaca       Ldb3
    Abhd3   Ugdh    Cbr2
    Cirbp   Iglon5  Lipe
    Mast1   C1qc    Hnrnpa3
    Slc44a2 Eif2b4  Septin9
    Magi2   Ythdf3  Hnrnpf
    Spag7   Rps15   Hnrnpab
    C1qbp   Atg4b   Sfxn3
    Pllp    Pvalb   Ctsz
    Smad3   Lig1    Coro1c
    Exoc6b  Ehd1    S100a11
    Chtop   Fn3krp  Ssrp1
    Hmgn2   Cxadr   Ptbp1
    Abhd6   Tmem186 Rac2
    Cers1   Lsm4    Nup85
    Dsg2    Rpl10   Nherf1
    Arl8b   Ca1 Coro1a
    Rbck1   Mindy2  Atp6v1a
    Abhd4   Rcc1    Arpc1b
    Frrs1l  Sart1   Krt8
    Rnf214  H3-5    Ruvbl1
    Cdc42ep1    H1-0    Pck1
    Sf3a3   Clca3a1 Hnrnpc
    Metap2  Lpcat4  Rida
    Eif2b4  H2-K1   Ces1d
    Inpp5f  S100a13 Smarcc2
    Urod    Trappc4 Ewsr1
    Nmt2    Ca4 Akr1c21
    Alg11   S100a4  Serpina3m
    Rad23b  Hbb-b2  Chmp1a
    Drap1   Celf1   Actbl2
    Hddc2   Rap2a   Ndufa1
    App Ca8 Ndufs4
    Cpsf6   Fgg Rala
    Ppp2r5c Ube2fb  C8a
    Kiss1   Eif3k   Mug1
    Rgs6    Ctnnbl1 Rrp12
    Txndc9  Timm10  mt-Nd2
    Ehd4    Uvrag   Apoa2
    Slc25a31    Arid1a  Ndufaf2
    Cpped1  Krt8    Mtx1
    Rbm27   Bpifb9a N6amt1
    Cmtm5   Col1a2  Mrpl49
    Gh1 Actbl2  Ndufs6
    Elp3    Evl Cisd1
    Atp6v0c Rpl39   Ndufa12
    Ints3   Col6a1  Serpina3k
    Tmem68  Map4k1
    Tiam1   Atp6v1e2
    Cox7a2l Col1a1
    Arfgef3 Rpl35
    Higd1a  Cpsf6
    Ptgds   Lum
    Eny2
    Tmsb10
    Per2
    Tatdn1
    Rpl30
    Irag2
    Otud6b
    Myl1
    Akap10
    Cav1
    Fbxl20
    Adh1
    Ehhadh
    Atp4a
    Naga
    Tceal3
    Col1a1
    Celf1
    Cops6
    Vps11
    Dync1i2
    Hm13
    Eif2s2
    Fgg
    Ralgapa1
    Ipo7
    Tmod3
    Fgb
    Phpt1
    Eif3j2
    Lars2
    Snrpb2
    Mug1
    Slc4a1
    Arg1
    Celf2
    Actbl2
    Krt6a
    Fga
    Ufd1
    Fam136a
    Ntmt1
    Opalin
    Omp
    Acsf3
    Lig3
    Xpo5
    Tmem240
    Aldob
    Adprs
    Stx1a
    H2bc14
    Serpina1d
    Cldn11
    Cyp2j5
    Slc6a11
    Cpm
    Acp2
    Ca1
    Cd44
    Cplx2")
fld <- fromList(as.list(d))
upset(fld)

As shown in the figure below, only five tissues could be seen, with the exception of the kidney, even in the presence of overlapping genes between the kidney and other tissues.
.
Can you kindly assist in this regard, I do not know if there is something wrong with my input data or the script itself.

No idea how your code looks like before pasting into the web form, but in general it is not a great idea pasting +260 rows/lines of data in your script.
From what I can see in the form you have columns of data separated by different number of spaces even in the header (2 and 4 spaces).
Because of this it is impossible to figure out to which column is which.
The easiest is to keep the data in the TSV or CSV format file and make sure there are no issues with column delimiters.
EDIT
I have saved the Google Spreadsheet as de_genes.csv, edited just in case Olfactory Bulb to Olfactory_Bulb then did:
in_fn <- "de_genes.csv"

df = fread(in_fn)

de_lists <- fromList(as.list(df))

upset(de_lists, sets=c( "Brainstem",    "Cerebellum",   "Olfactory_Bulb",   "Heart", "Kidney", "Liver") , order.by = "freq")

I got:

I tried using that strategy (i.e., keeping data in csv instead of pasting it) but I am struggling with reading my csv file using the "data.table" function.        See the original data here (https://docs.google.com/spreadsheets/d/17KuGuzBjBhaRBoDCq-3SvnklBZQF1C-HcTpss1V0zKA/edit?usp=sharing)

Updated the answer. It should work now

It is indeed working, thank you.

@S'boneloGlenKhumato  could you upvote this answer and mark it as accepted please?


text:
Multi Factor in Deseq2 Gene enrichment analysis
I want to see how the gene expression differs in breast cancer between three species, and I am using DESeqDataSetFromMatrix on my count table.
dds <- DESeqDataSetFromMatrix(countData = whole_data, colData = whole_col_data, design = ~ species)
dds$species <- factor(dds$species, levels=c("human","canine","mouse"))
dds$species <- relevel(dds$species, ref = "human")
dds <- DESeq(dds)

Although I specified my levels for all three species my results are still as follows and only compare between two species:
log2 fold change (MLE): species **mouse vs human**
Wald test p-value: species **mouse vs human**
DataFrame with 13556 rows and 6 columns
        baseMean log2FoldChange     lfcSE      stat      pvalue

Can someone please explain to me why that is happening?

If you don't tell DESeq what contrast to make, it just does the default contrast.  Which is your reference versus the last other one alphabetically.

Oh understood, Thanks! I thought I was doing that with the factor command, can you please let me know ow should I do that, is there a tutorial you recommend for it?

You need results(dds, contrast = c("species", "human", "canine")) if you're looking to compare human vs canine. Three params, first is the variable (column) name, 2 and 3 are the conditions you're comparing. See: https://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#note-on-factor-levels

Thank you Ram RS.

I just noticed that in contrast there are juts humans and canines still, I want the three species to be there, human, canine, and mouse. Since I have three different species and want to compare the expression among them. I added the mouse to the contrast but I got an error "contrast', as a character vector of length 3, should have the form: contrast = c('factorName','numeratorLevel','denominatorLevel'),
see the manual page of ?results for more information" is it even possible to do this for more than two categories?

You can only compare two conditions. You can of course run the results multiple times with different pairs of conditions but you cannot compare 3 groups using a 2-group statistical test.

Thanks, Ram RS!

Wouldn't it be better to run the DESeqDataSetFromMatrix on a matrix that has only samples of two conditions then? instead of giving all the categories in the same input matrix? I mean I can give the count table of mouse and human, human and canine, and canine and mouse separately if it's not considering them all, would that make a change? Thanks for helping!

In general, it is better not to subset the data.  Keep it all together, use contrasts to specify what to compare to what.

Thanks for your answer! Just out of curiosity do you know the main reason behind it?

Better dispersion estimates. Part of the issue was the RNA-Seq is that people generally don't have the money to do all the replicates that are really required for robust results.  Part of what DESeq does is to squeeze as much statistical information out as possible of the data you do have.  Removing real sample data hurts that effort, usually.

source: Stack Exchange

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

sort: votes


text:
Can alphafold 2 be used to identify residues important for binding?
I am interested in using AlphaFold2 to identify residues within a homodimer that are important for the protein-protein interaction. I have been able to use ColabFold to produce a predicted dimer structure - however, as far as I can tell AlphaFold does not actually identify the residues or regions that are actually important for binding. Is there a way to manipulate the software to do this?
Thanks!

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


text:
Assertion `compressBound(0xff00) < 0x10000' failed when using STAR 2.7.9a
This question was also asked on github
I am using STAR 2.7.9a. In the last two days, I am getting an error when running alignments with STAR (independently from the data and genome index I am using). What could be causing this error, which is reported in the next paragraph? Notably, STAR returns this error even when performing alignments that used to complete successfully in the past, before this issue started. Thank you very much for your attention and help on this!
INFO: localdir at /localscratch/8889241.roberto.albanese
INFO: your job will run on local system.
STAR: bgzf.c:158: bgzf_open: Assertion `compressBound(0xff00) < 0x10000' failed.
*** Aborted
Register dump:

 RAX: 0000000000000000   RBX: 0000000000000006   RCX: 00007ffff684a37f
 RDX: 0000000000000000   RSI: 00007ffffffdb580   RDI: 0000000000000002
 RBP: 00007ffff699d5f8   R8 : 0000000000000000   R9 : 00007ffffffdb580
 R10: 0000000000000008   R11: 0000000000000246   R12: 0000000000594eab
 R13: 0000000000594f18   R14: 000000000000009e   R15: 00007ffffffeb540
 RSP: 00007ffffffdb580

 RIP: 00007ffff684a37f   EFLAGS: 00000246

 CS: 0033   FS: 0000   GS: 0000

 Trap: 00000000   Error: 00000000   OldMask: 00000000   CR2: 00000000

 FPUCW: 0000037f   FPUSW: 00000000   TAG: 00007fff
 RIP: f76a5f3b   RDP: 00000000

 ST(0) 0000 0000000000000000   ST(1) 0000 0000000000000000
 ST(2) 0000 0000000000000000   ST(3) 0000 0000000000000000
 ST(4) 0000 0000000000000000   ST(5) 0000 0000000000000000
 ST(6) ffff b000000000000000   ST(7) b000 b000000000000000
 mxcsr: 1fa0
 XMM0:  00000000000000000000000000000000 XMM1:  00000000000000000000000000000000
 XMM2:  00000000000000000000000000000000 XMM3:  00000000000000000000000000000000
 XMM4:  00000000000000000000000000000000 XMM5:  00000000000000000000000000000000
 XMM6:  00000000000000000000000000000000 XMM7:  00000000000000000000000000000000
 XMM8:  00000000000000000000000000000000 XMM9:  00000000000000000000000000000000
 XMM10: 00000000000000000000000000000000 XMM11: 00000000000000000000000000000000
 XMM12: 00000000000000000000000000000000 XMM13: 00000000000000000000000000000000
 XMM14: 00000000000000000000000000000000 XMM15: 00000000000000000000000000000000

Backtrace:
/lib64/libc.so.6(gsignal+0x10f)[0x7ffff684a37f]
/lib64/libc.so.6(abort+0x127)[0x7ffff6834db5]
/lib64/libc.so.6(+0x21c89)[0x7ffff6834c89]
/lib64/libc.so.6(+0x2fa76)[0x7ffff6842a76]
/home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR[0x53b3b0]
/home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR[0x4af101]
/home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR[0x40da52]
/lib64/libc.so.6(__libc_start_main+0xf3)[0x7ffff6836493]
/home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR[0x419125]

Memory map:

00400000-005bc000 r-xp 00000000 00:48 4804921705                         /home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR
007bc000-007bd000 r--p 001bc000 00:48 4804921705                         /home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR
007bd000-007cb000 rw-p 001bd000 00:48 4804921705                         /home/roberto.albanese/software/STAR-2.7.9a/bin/Linux_x86_64/STAR
007cb000-00832000 rw-p 00000000 00:00 0                                  [heap]
7ffff660f000-7ffff6612000 r-xp 00000000 08:02 36793927                   /usr/lib64/libdl-2.28.so
7ffff6612000-7ffff6811000 ---p 00003000 08:02 36793927                   /usr/lib64/libdl-2.28.so
7ffff6811000-7ffff6812000 r--p 00002000 08:02 36793927                   /usr/lib64/libdl-2.28.so
7ffff6812000-7ffff6813000 rw-p 00003000 08:02 36793927                   /usr/lib64/libdl-2.28.so
7ffff6813000-7ffff69cf000 r-xp 00000000 08:02 36776071                   /usr/lib64/libc-2.28.so
7ffff69cf000-7ffff6bce000 ---p 001bc000 08:02 36776071                   /usr/lib64/libc-2.28.so
7ffff6bce000-7ffff6bd2000 r--p 001bb000 08:02 36776071                   /usr/lib64/libc-2.28.so
7ffff6bd2000-7ffff6bd4000 rw-p 001bf000 08:02 36776071                   /usr/lib64/libc-2.28.so
7ffff6bd4000-7ffff6bd8000 rw-p 00000000 00:00 0
7ffff6bd8000-7ffff6bf3000 r-xp 00000000 08:02 37169965                   /usr/lib64/libpthread-2.28.so
7ffff6bf3000-7ffff6df2000 ---p 0001b000 08:02 37169965                   /usr/lib64/libpthread-2.28.so
7ffff6df2000-7ffff6df3000 r--p 0001a000 08:02 37169965                   /usr/lib64/libpthread-2.28.so
7ffff6df3000-7ffff6df4000 rw-p 0001b000 08:02 37169965                   /usr/lib64/libpthread-2.28.so
7ffff6df4000-7ffff6df8000 rw-p 00000000 00:00 0
7ffff6df8000-7ffff6e0f000 r-xp 00000000 08:02 100684981                  /cm/local/apps/gcc/10.2.0/lib64/libgcc_s.so.1
7ffff6e0f000-7ffff700e000 ---p 00017000 08:02 100684981                  /cm/local/apps/gcc/10.2.0/lib64/libgcc_s.so.1
7ffff700e000-7ffff700f000 r--p 00016000 08:02 100684981                  /cm/local/apps/gcc/10.2.0/lib64/libgcc_s.so.1
7ffff700f000-7ffff7010000 rw-p 00017000 08:02 100684981                  /cm/local/apps/gcc/10.2.0/lib64/libgcc_s.so.1
7ffff7010000-7ffff704e000 r-xp 00000000 08:02 100684988                  /cm/local/apps/gcc/10.2.0/lib64/libgomp.so.1.0.0
7ffff704e000-7ffff724d000 ---p 0003e000 08:02 100684988                  /cm/local/apps/gcc/10.2.0/lib64/libgomp.so.1.0.0
7ffff724d000-7ffff724e000 r--p 0003d000 08:02 100684988                  /cm/local/apps/gcc/10.2.0/lib64/libgomp.so.1.0.0
7ffff724e000-7ffff724f000 rw-p 0003e000 08:02 100684988                  /cm/local/apps/gcc/10.2.0/lib64/libgomp.so.1.0.0
7ffff724f000-7ffff73d0000 r-xp 00000000 08:02 36888910                   /usr/lib64/libm-2.28.so
7ffff73d0000-7ffff75cf000 ---p 00181000 08:02 36888910                   /usr/lib64/libm-2.28.so
7ffff75cf000-7ffff75d0000 r--p 00180000 08:02 36888910                   /usr/lib64/libm-2.28.so
7ffff75d0000-7ffff75d1000 rw-p 00181000 08:02 36888910                   /usr/lib64/libm-2.28.so
7ffff75d1000-7ffff7794000 r-xp 00000000 08:02 100685649                  /cm/local/apps/gcc/10.2.0/lib64/libstdc++.so.6.0.28
7ffff7794000-7ffff7993000 ---p 001c3000 08:02 100685649                  /cm/local/apps/gcc/10.2.0/lib64/libstdc++.so.6.0.28
7ffff7993000-7ffff799e000 r--p 001c2000 08:02 100685649                  /cm/local/apps/gcc/10.2.0/lib64/libstdc++.so.6.0.28
7ffff799e000-7ffff79a1000 rw-p 001cd000 08:02 100685649                  /cm/local/apps/gcc/10.2.0/lib64/libstdc++.so.6.0.28
7ffff79a1000-7ffff79a4000 rw-p 00000000 00:00 0
7ffff79a4000-7ffff79ca000 r-xp 00000000 00:89 7817491573                 /share/sciswrepo/spack/latest-2024-01-17/linux-centos8-skylake_avx512/gcc-8.4.1/zlib-ng-2.0.0-qvhoazdgq5hmo6s32psbvmz75vofgtbt/lib/libz.so.1.2.11.zlib-ng
7ffff79ca000-7ffff7bc9000 ---p 00026000 00:89 7817491573                 /share/sciswrepo/spack/latest-2024-01-17/linux-centos8-skylake_avx512/gcc-8.4.1/zlib-ng-2.0.0-qvhoazdgq5hmo6s32psbvmz75vofgtbt/lib/libz.so.1.2.11.zlib-ng
7ffff7bc9000-7ffff7bca000 r--p 00025000 00:89 7817491573                 /share/sciswrepo/spack/latest-2024-01-17/linux-centos8-skylake_avx512/gcc-8.4.1/zlib-ng-2.0.0-qvhoazdgq5hmo6s32psbvmz75vofgtbt/lib/libz.so.1.2.11.zlib-ng
7ffff7bca000-7ffff7bcb000 rw-p 00026000 00:89 7817491573                 /share/sciswrepo/spack/latest-2024-01-17/linux-centos8-skylake_avx512/gcc-8.4.1/zlib-ng-2.0.0-qvhoazdgq5hmo6s32psbvmz75vofgtbt/lib/libz.so.1.2.11.zlib-ng
7ffff7bcb000-7ffff7bcf000 r-xp 00000000 08:02 36760911                   /usr/lib64/libSegFault.so
7ffff7bcf000-7ffff7dce000 ---p 00004000 08:02 36760911                   /usr/lib64/libSegFault.so
7ffff7dce000-7ffff7dcf000 r--p 00003000 08:02 36760911                   /usr/lib64/libSegFault.so
7ffff7dcf000-7ffff7dd0000 rw-p 00004000 08:02 36760911                   /usr/lib64/libSegFault.so
7ffff7dd0000-7ffff7dfc000 r-xp 00000000 08:02 36760102                   /usr/lib64/ld-2.28.so
7ffff7fd7000-7ffff7fdd000 rw-p 00000000 00:00 0
7ffff7fef000-7ffff7ff6000 rw-p 00000000 00:00 0
7ffff7ff6000-7ffff7ffa000 r--p 00000000 00:00 0                          [vvar]
7ffff7ffa000-7ffff7ffc000 r-xp 00000000 00:00 0                          [vdso]
7ffff7ffc000-7ffff7ffd000 r--p 0002c000 08:02 36760102                   /usr/lib64/ld-2.28.so
7ffff7ffd000-7ffff7fff000 rw-p 0002d000 08:02 36760102                   /usr/lib64/ld-2.28.so
7ffffffcd000-7ffffffff000 rw-p 00000000 00:00 0                          [stack]
ffffffffff600000-ffffffffff601000 r-xp 00000000 00:00 0                  [vsyscall]
/cm/local/apps/slurm/var/spool/job8889241/slurm_script: line 18: 3742255 Aborted                 (core dumped)
```

bgzf.c:158: bgzf_open: Assertion `compressBound(0xff00) < 0x10000' failed.

Upstream HTSlib users encountered this problem when HTSlib was used in conjunction with zlib-ng, which has slightly worse worst-case compression behaviour than the original zlib.
Is zlib-ng installed on your machine, and is STAR using it rather than vanilla zlib? (If so, presumably when these alignments completed successfully in the past, that was before zlib-ng was installed.)
Whether the problem is due to zlib-ng or you have managed in some other way to have a slightly too-big compressed block for some particular input data, the same zlib-ng-supporting fix will resolve the issue because it removed these assertions and added code to handle any too-big blocks encountered.
This problem was fixed in HTSlib three years ago.
Unfortunately STAR bundles a pre-historic copy of HTSlib that predates this fix. You have also reported this as a STAR issue, so hopefully the authors will address it in due course.
In the meantime, you could for example use bioconda's STAR package, which is built to use bioconda's HTSlib package, which is up-to-date and does not exhibit this assertion failure. Or, if you are aware of having installed zlib-ng about two days ago, you could uninstall it.


text:
Aside from the header info, is there any way to distinguish between sequences generated by MiSeq and MiniSeq platforms?
I am comparing Illumina MiSeq (using MiSeq v3 kit) and MiniSeq (MiniSeq High Output kit) instruments for sequencing E.coli and Salmonella. I am curious if there is any way to bioinformatically distinguish between the sequences generated by these instruments without the header information? Would it simply be the read length difference (MiSeq=2 x 300bp vs MiniSeq=2 x 150 bp) and output (MiSeq=13.2-15Gb vs MiniSeq=6.6-7.5Gb). Quality wise they look nearly identical in my experience.
Thank you.

Can you explain why? I mean, why would you need this? Presumably, you know what sequencer you used, right? And if not, you can use the header. Why would you need to hobble yourself this way? And why aren't the differences you listed (one produces reads of twice the length and files of twice the size of the other) already enough? What else are you looking for?

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Why could you possibly need this?? MiSeq is higher throughput so you should have more reads under certain assumptions


text:
How to parse hmmsearch output?
I did hmmsearch with --notextw and --domtblout and I see several results for ids, I want to get the best result for each id. How I can parse it?

I saw the tools in Biopython SearchIO but I really can't understand how I can do it. Please, can you explain it in detail?
This is part of my file:
NODE_1_length_959108_cov_67.987603_423 -            531 2-Hacid_dh           PF00389.34   134   1.6e-35  120.1   0.4   1   1   1.1e-38   3.3e-35  119.2   0.4     1   134     9   316     9   316 0.99 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_7_length_349797_cov_70.130335_175 -            319 2-Hacid_dh           PF00389.34   134   1.3e-13   49.2   0.0   1   1   5.1e-17   1.5e-13   49.0   0.0    21   129    28   314     6   318 0.94 # 187752 # 188708 # -1 # ID=7_175;partial=00;start_type=GTG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.637
NODE_1_length_959108_cov_67.987603_423 -            531 2-Hacid_dh_C         PF02826.23   178   6.6e-63  209.7   0.0   1   2   1.3e-65   1.9e-62  208.2   0.0     1   178   110   284   110   284 0.98 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_1_length_959108_cov_67.987603_423 -            531 2-Hacid_dh_C         PF02826.23   178   6.6e-63  209.7   0.0   2   2       1.1   1.6e+03   -3.3   0.0     6    44   331   369   328   371 0.66 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_7_length_349797_cov_70.130335_175 -            319 2-Hacid_dh_C         PF02826.23   178   7.2e-51  170.5   0.0   1   1   6.8e-54     1e-50  170.0   0.0     7   178   118   287   112   287 0.94 # 187752 # 188708 # -1 # ID=7_175;partial=00;start_type=GTG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.637
NODE_7_length_349797_cov_70.130335_232 -            309 2-Hacid_dh_C         PF02826.23   178   2.2e-49  165.7   0.0   1   1   1.9e-52   2.9e-49  165.3   0.0     4   178    97   271    95   271 0.92 # 243113 # 244039 # -1 # ID=7_232;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.653
NODE_9_length_318862_cov_66.931761_127 -            311 2-Hacid_dh_C         PF02826.23   178   3.8e-48  161.7   0.0   1   1   3.3e-51   4.9e-48  161.3   0.0     2   176   100   271    98   273 0.89 # 151217 # 152149 # -1 # ID=9_127;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.652
NODE_2_length_853960_cov_68.957247_395 -            506 2-oxoacid_dh         PF00198.27   232   1.7e-82  274.8   0.0   1   1   2.2e-85   2.3e-82  274.4   0.0     2   231   277   504   276   505 0.98 # 406642 # 408159 # 1 # ID=2_395;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=5-10bp;gc_cont=0.644
NODE_1_length_959108_cov_67.987603_423 -            531 PGDH_inter           PF19304.3    119   9.9e-19   66.1   1.5   1   3      0.76   4.6e+03   -3.8   0.0    68    81    87   100    86   123 0.54 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_1_length_959108_cov_67.987603_423 -            531 PGDH_inter           PF19304.3    119   9.9e-19   66.1   1.5   2   3       0.4   2.4e+03   -2.9   0.0    54    72   208   226   199   249 0.73 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_1_length_959108_cov_67.987603_423 -            531 PGDH_inter           PF19304.3    119   9.9e-19   66.1   1.5   3   3   2.6e-22   1.6e-18   65.4   0.4     2   119   326   446   325   446 0.91 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_2_length_853960_cov_68.957247_550 -            554 PGI                  PF00342.23   486  8.5e-202  669.5   0.0   1   1  1.6e-205  9.7e-202  669.3   0.0     1   486    52   543    52   543 0.97 # 568442 # 570103 # -1 # ID=2_550;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.626

You can see NODE_1_length_959108_cov_67.987603_423 several times with different annotations (2-Hacid_dh, PF00389.34 and 2-Hacid_dh_C, PF02826.23 and PGDH_inter, PF19304.3).
I need to get only best result, only ONE result. I need to see only one row with this id (NODE_1_length_959108_cov_67.987603_423) with  only one annotation.
I assumed that there were some programs with certain algorithms that could select one thing. For example, as in blastp there is --max_hsps 1. I need to make a pivot table for each ID. I did a search in the Swissprot database, in the Pataa database, and also in Pfam. But there are several search results in pfam and I don’t know how to put it in the pivot table

Can you please paste the content as text instead of using screenshots? Also, the user manual says something really weird: "The tblout format is deliberately space-delimited (rather than tab-delimited) and justified into aligned columns, so these files are suitable both for automated parsing and for human examination." I think the author is wrong but your output does not look justified either, so I don't know what's happening.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Please post the real file. We need an example we can use. Your image shows headers, but they are not shown in the text example. Please make sure the example accurately reflects the real input, and then also add what output you expect. You say you want the "best" hit. OK, but "best" based on what? Lowest E-value? Highest score? Longest length? And why would you only want one row? What if many are equally good?

blast's --max_hsps 1 does not work the way you think it does, it does not return the best hit.

Please [edit] and give us these details. If we understand the biological question you are trying to answer, we might be able to help more. We don't even know the species you are working with or what input you used. You mentioned Pfam, should we assume you are looking for protein domain matches? To find homologs in different species perhaps? If so, please tell us since that changes what "best" would mean. Do you want the best match to the search HMM? The longest match for each input? And please add a full example. We need to see the headers in order to parse and to know what the data are.

Hi @bobr could kindly mark one of the answers as "accepted" please? The question has been well answered - so its only fair on those answering - in addition it helps site stats (the question can keep circulating within the site otherwise).

The TSV output should be sorted by default in a way that for a given query the top hit comes first, followed by other hits. For a quick extraction of top hits you may use a rather simplistic script:
https://github.com/Cancer-Genomics-TH/tRNA_methyl/blob/main/src/parse_top_hits_cmscan.py

Since we don't know what you are trying to maximize here, which of the various columns would define the "best" hit for what you are trying to do, here is a generic solution that will just find the highest score:
$ awk '( ! best[$1] ||  $6 > best[$1] ){ best[$1]=$6; line[$1]=$0 } END{ for(name in best){ print line[name] }}' file
NODE_1_length_959108_cov_67.987603_423 -            531 2-Hacid_dh_C         PF02826.23   178   6.6e-63  209.7   0.0   1   2   1.3e-65   1.9e-62  208.2   0.0     1   178   110   284   110   284 0.98 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_7_length_349797_cov_70.130335_232 -            309 2-Hacid_dh_C         PF02826.23   178   2.2e-49  165.7   0.0   1   1   1.9e-52   2.9e-49  165.3   0.0     4   178    97   271    95   271 0.92 # 243113 # 244039 # -1 # ID=7_232;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.653
NODE_7_length_349797_cov_70.130335_175 -            319 2-Hacid_dh_C         PF02826.23   178   7.2e-51  170.5   0.0   1   1   6.8e-54     1e-50  170.0   0.0     7   178   118   287   112   287 0.94 # 187752 # 188708 # -1 # ID=7_175;partial=00;start_type=GTG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.637
NODE_9_length_318862_cov_66.931761_127 -            311 2-Hacid_dh_C         PF02826.23   178   3.8e-48  161.7   0.0   1   1   3.3e-51   4.9e-48  161.3   0.0     2   176   100   271    98   273 0.89 # 151217 # 152149 # -1 # ID=9_127;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.652
NODE_2_length_853960_cov_68.957247_395 -            506 2-oxoacid_dh         PF00198.27   232   1.7e-82  274.8   0.0   1   1   2.2e-85   2.3e-82  274.4   0.0     2   231   277   504   276   505 0.98 # 406642 # 408159 # 1 # ID=2_395;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=5-10bp;gc_cont=0.644
NODE_2_length_853960_cov_68.957247_550 -            554 PGI                  PF00342.23   486  8.5e-202  669.5   0.0   1   1  1.6e-205  9.7e-202  669.3   0.0     1   486    52   543    52   543 0.97 # 568442 # 570103 # -1 # ID=2_550;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.626

To get the maximum value of a different field, just change $6 to the corresponding field number.
And here's a version that looks for the smallest value, in this case using the E-value (field 5):
$ awk '( ! best[$1] ||  $6 < best[$1] ){ best[$1]=$6; line[$1]=$0 } END{ for(name in best){ print line[name] }}' file
NODE_1_length_959108_cov_67.987603_423 -            531 PGDH_inter           PF19304.3    119   9.9e-19   66.1   1.5   1   3      0.76   4.6e+03   -3.8   0.0    68    81    87   100    86   123 0.54 # 448882 # 450474 # -1 # ID=1_423;partial=00;start_type=GTG;rbs_motif=GGAG/GAGG;rbs_spacer=5-10bp;gc_cont=0.664
NODE_7_length_349797_cov_70.130335_232 -            309 2-Hacid_dh_C         PF02826.23   178   2.2e-49  165.7   0.0   1   1   1.9e-52   2.9e-49  165.3   0.0     4   178    97   271    95   271 0.92 # 243113 # 244039 # -1 # ID=7_232;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.653
NODE_7_length_349797_cov_70.130335_175 -            319 2-Hacid_dh           PF00389.34   134   1.3e-13   49.2   0.0   1   1   5.1e-17   1.5e-13   49.0   0.0    21   129    28   314     6   318 0.94 # 187752 # 188708 # -1 # ID=7_175;partial=00;start_type=GTG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.637
NODE_9_length_318862_cov_66.931761_127 -            311 2-Hacid_dh_C         PF02826.23   178   3.8e-48  161.7   0.0   1   1   3.3e-51   4.9e-48  161.3   0.0     2   176   100   271    98   273 0.89 # 151217 # 152149 # -1 # ID=9_127;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.652
NODE_2_length_853960_cov_68.957247_395 -            506 2-oxoacid_dh         PF00198.27   232   1.7e-82  274.8   0.0   1   1   2.2e-85   2.3e-82  274.4   0.0     2   231   277   504   276   505 0.98 # 406642 # 408159 # 1 # ID=2_395;partial=00;start_type=ATG;rbs_motif=GGA/GAG/AGG;rbs_spacer=5-10bp;gc_cont=0.644
NODE_2_length_853960_cov_68.957247_550 -            554 PGI                  PF00342.23   486  8.5e-202  669.5   0.0   1   1  1.6e-205  9.7e-202  669.3   0.0     1   486    52   543    52   543 0.97 # 568442 # 570103 # -1 # ID=2_550;partial=00;start_type=ATG;rbs_motif=None;rbs_spacer=None;gc_cont=0.626

And here is the same code, expanded, with an explanation (this will still execute correctly if you copy/paste it into a terminal):
awk '
  ## If there's no value for the first field in the 'best' array or, if the
  ## existing value is lower than the 6th field, then save the $6th  field
  ## as the new value for this name in the 'best' array, and save the current
  ## line as the value in the 'line' array.
  ( ! best[$1] ||  $6 < best[$1] ){
     best[$1]=$6
     line[$1]=$0
  }
  ## When we reach the end of the file, iterate over the 'line'
  ## array, and print each element.
  END{
    for(name in best){
      print line[name]
    }
   }' file


text:
Single Node Deletion in Cheng and Church Algorithm
Cheng and Church Algorithm is a biclustering algorithm or an unsupervised learning task with which both rows and columns in a data matrix can be placed in concurrent clusters.
In the original paper, authors explained that a row or column of an input matrix $A$ with the largest $d(i)$ or $d(j)$ should be removed in single node deletion phase:

Find the row $i \in I$ with largest: $d(i) = \frac{\sum_{j \in J} (a_{ij} - a_{iJ} - a_{Ij} + a_{IJ})}{\lvert J \rvert}$
Find the column $j \in J$ with largest: $d(j) = \frac{\sum{i \in I} (a_{ij} - a_{iJ} - a_{Ij} + a_{IJ})}{\lvert I \rvert}$

I found one implementation of the algorithm trying to do the following operations for single node deletion:
row_max_msr = np.argmax(row_msr) #finding the index of rows having largest d(i)
col_max_msr = np.argmax(col_msr) #finding the index of columns having largest d(j)

#Continue with removing row where $d(i) >= d(j)$ or column $d(j) > d(i)$
if row_msr[row_max_msr] >= col_msr[col_max_msr]:
    row2remove = row_indices[row_max_msr]
    rows[row2remove] = False
else:
    col2remove = col_indices[col_max_msr]
    cols[col2remove] = False

I was wondering whether anyone has a suggestion for this matter. I don't mean that the code has installation or execution problem. Actually, I want to be sure of the logic to be aligned with the algorithm.

Padilha, V. A. & Campello, R. J. G. B. (2017). A systematic comparative evaluation of biclustering techniques. BMC Bioinformatics, 18(1):55.

https://github.com/padilha/biclustlib
It needs the r-bicluster library which is available on conda https://anaconda.org/conda-forge/r-biclust
Citing the authors GitHub page for instructions
import numpy as np

from biclustlib.algorithms import ChengChurchAlgorithm
from biclustlib.datasets import load_yeast_tavazoie

# load yeast data used in the original Cheng and Church's paper
data = load_yeast_tavazoie().values

# missing value imputation suggested by Cheng and Church
missing = np.where(data < 0.0)
data[missing] = np.random.randint(low=0, high=800, size=len(missing[0]))

# creating an instance of the ChengChurchAlgorithm class and running with the parameters of the original study
cca = ChengChurchAlgorithm(num_biclusters=100, msr_threshold=300.0, multiple_node_deletion_threshold=1.2)
biclustering = cca.run(data)
print(biclustering)

Thank you for the response, however, I don't mean that the code has installation or execution problem. Actually, I want to be sure of the logic to be aligned with the algorithm.


text:
Calculating Fisher's exact test for COG categories
This is a continuation of the following question -
Fisher's exact test for COG categories in pan, core genome analysis
Also, the dataset -

How do I perform a Fisher's exact test on this data to show that a COG category is significantly more or less in the pan/core/accessory/species specific genomes?

I'll post a response late Monday, sorry really busy. Basically you can do spot checks via Fishers Exact test or go the whole hog and do a Chi2 test. The matrix you have is pan/core vs cogs.

Hi @M__, good to know. Thanks a lot!

What I would personally do is run a Χ2
from scipy.stats import chisquare
import numpy as np

pan = [11, 39, 4 ....]
core = [29, 50, 7 ... ]
f_exp = np.array(pan)
f_obs = np. array(core)
print (chisquare(f_obs=f_obs, f_exp=f_exp))

This will output the total statistic and whether its above or below the critical probability.
Fisher's Exact test is essentially Χ2 for a small number of observations, which is not the case with your data.
Looking at your data Χ2 is not likely to be "significant" (which is the result you are looking for), i.e. the pan genome is the core genome across all COGs.

Thank you very much for all your efforts @M__ ! Hopefully, this works.


text:
How to debug a bash script
Please somebody debug my code
In sub-folders there are 500000 vcf files from which I only need 500 vcf files listed in temp.txt which randomly distributed in sub-folders 10 to 60. The code reads each of 500000 vcf files one by one, if a file was in temp.txt, the code runs
dx run swiss-army-knife -icmd="tabix -h $vcf -R $REGIONS_FILE > $output_file" -iin="$filename_root" -iin="$vcf.*" -iin="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/file.bed" --destination "$OUTPUT_FOLDER" --detach --priority high --yes
    done < "$VCF_LIST_FILE"

on that if not, the code goes to next file such a process takes years. I want this code only locates files listed in temp.txt and run the main code on them by iteration.
#!/bin/bash

# Set the working directory to the desired location within DNAnexus

# Number of concurrent processes
N=4

# Iterate from 10 to 60
for folder_number in {10..60}; do
  (
    # Define the input files and parameters for this iteration
    folder_path="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/$folder_number"
    VCF_LIST_FILE="temp.txt"
    REGIONS_FILE='file.bed'
    OUTPUT_FOLDER='/out'  # Replace with your desired output folder

    # Set the working directory to the folder for this iteration
    dx cd "$folder_path"

    # Read each VCF file from the list
    while read -r vcf; do
      filename_root=$(basename "$vcf")
      output_file="${filename_root}_extracted.vcf.gz"

      # Use dx run swiss-army-knife to execute tabix on the VCF file
      dx run swiss-army-knife -icmd="tabix -h $vcf -R $REGIONS_FILE > $output_file" -iin="$filename_root" -iin="$vcf.*" -iin="/Bulk/DRAGEN\ WGS/Whole\ genome\ variant\ call\ files\ (VCFs)\ (DRAGEN)\ [500k release]/file.bed" --destination "$OUTPUT_FOLDER" --detach --priority high --yes
    done < "$VCF_LIST_FILE"
  ) &

  # Limit the number of concurrent processes
  ((i=i%N)); ((i++==0)) && wait
done

# Wait for all background jobs to finish
wait

# Finished message
echo "++ Finished extracting genomic coordinates from VCF files"

This is how temp.txt looks. The first two digits of vcf file names defines which sub-folder they are located. For instance, the below vcf files distributed in sub-folders 32, 33 and 34
head temp_new.txt
head: temp_new.txt: No such file or directory
izadia@LV19Y7325V Downloads % head new_temp.txt
11/1195100_24053_0_0.dragen.hard-filtered.vcf.gz
12/1204363_24053_0_0.dragen.hard-filtered.vcf.gz
12/1246629_24053_0_0.dragen.hard-filtered.vcf.gz
13/1311312_24053_0_0.dragen.hard-filtered.vcf.gz
14/1427973_24053_0_0.dragen.hard-filtered.vcf.gz
14/1435043_24053_0_0.dragen.hard-filtered.vcf.gz
15/1558824_24053_0_0.dragen.hard-filtered.vcf.gz
16/1683316_24053_0_0.dragen.hard-filtered.vcf.gz
16/1692707_24053_0_0.dragen.hard-filtered.vcf.gz
17/1749825_24053_0_0.dragen.hard-filtered.vcf.gz
izadia@LV19Y7325V Downloads % tail new_temp.txt
52/5230856_24053_0_0.dragen.hard-filtered.vcf.gz
52/5257892_24053_0_0.dragen.hard-filtered.vcf.gz
54/5453459_24053_0_0.dragen.hard-filtered.vcf.gz
55/5514349_24053_0_0.dragen.hard-filtered.vcf.gz
55/5526010_24053_0_0.dragen.hard-filtered.vcf.gz
57/5718723_24053_0_0.dragen.hard-filtered.vcf.gz
58/5868836_24053_0_0.dragen.hard-filtered.vcf.gz
58/5875882_24053_0_0.dragen.hard-filtered.vcf.gz
58/5889611_24053_0_0.dragen.hard-filtered.vcf.gz
59/5917345_24053_0_0.dragen.hard-filtered.vcf.gz
izadia@LV19Y7325V Downloads %

After running the script I get this error
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1195100_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1204363_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1246629_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1311312_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1427973_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1435043_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1558824_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1683316_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1692707_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID
dxpy.exceptions.DXCLIError: Value provided for input field "in" could not be parsed as file: could not resolve "1749825_24053_0_0.dragen.hard-filtered.vcf.gz" to a name or ID

because the codes starts searching from sub-folder 10 for each of files in new_temp.txt one by one

So what do you want to do? If you don't want to check if the files exist, then don't check if the files exist. Are there any vcf files in your various folders that should not be processed? What is the point of temp,txt? And what does it contain? Full paths? FIle names?

Oh, and your code never checks if a file exists. Just [edit] your question, explain what you really want to do, and maybe we can help.

I edited my main post

And what does temp.txt look like? Please show us an example. We need to know if it is one file name per line, many per line, if it is just the names or the full paths etc.

While you can do it in bash often switching to ie python gives you more debugging options, cleaner code, etc.

Please do not edit your posts as follow up to answers without making the edits extremely clear to follow. Right now, we need to read the answers and your edit history (not question) side by side to understand what's going on and that's a huge pain.

Please don't forget to mark the answer below as "accepted" if this addresses your question (it does look like it).

Thanks for the answer but the problem have not been solved yet so the code still stats checking from sub-folder 10 for the listed files one by one which takes forever

The first two digits of vcf file names defines which sub-folder they are located.

You could use this information to edit temp.txt so it contains the full path:
grep -v ".tbi$" temp.txt | sed 's/^\(.\)\(.\)/\1\2\/\1\2/' > new_temp.txt

I'm assuming here that you don't actually want to run dx run swiss-army-knife on the vcf indexes.
This will turn the example of temp.txt you gave into this:
32/3295094_24053_0_0.dragen.hard-filtered.vcf.gz
33/3309510_24053_0_0.dragen.hard-filtered.vcf.gz
34/3468518_24053_0_0.dragen.hard-filtered.vcf.gz
34/3471960_24053_0_0.dragen.hard-filtered.vcf.gz

You can then simply run a single for loop over that, instead of the nested one you are using now.
eg.
while read vcf
  do
    dx run swiss-army-knife [...]
  done <  new_temp.txt


text:
Comparing sequence similarity of two sets of peptide sequences?
What tools can we use to compare two sets of peptide sequences?

The biological question is I have some kmer sequences of some proteins and would like to know how similar are these kmer sequences compared with known MHC classs I and II peptides.
Therefore I need to compare two sets of peptides sequences.
The kmer and peptide sequences lengths range from 5-80aa.

Alternatively, I am also considering to use the full length of the proteins for comparison.

What kind of comparison? All against all? Pairs? Triplets? Are you expecting them all to be similar? Are you trying to group them into sets based on similarity? Blastp doesn't seem relevant at all here (and 11 is the default word size of blastn not blastp, but that's neither here nor there). Please [edit] your question and tell us what you are trying to do, what is the biological question you are trying to answer. We need to understand that in order to suggest approaches.

Hi @terdon, thanks for the suggestion and now edited the question.

Just a reminder to consider marking one of the questions as "answered". This is definitely helpful for all concerned.

Not sure what is the shortest peptide one can fish out with LAST out of small protein database (well, it starts with protein FASTA file), but you may check out LAST cookbook. This can use different PAM matrices so you have the proper similar amino acid matching. If that does not suit your needs there can be some more desperate route using string fuzzy matching i.e. with thefuzz python library and checking the best Levenshtein distance scores.

I would suggest focusing on short k-mers if you are using very small sequences. Simply quantify the similarity of the k-mer profile of each relevant pair of peptides. There are many tools for analyzing k-mer profiles, though many are targeted towards DNA so you should keep an eye out for that. Here is one python tool.
It is also not difficult to write a k-mer counting script if you decide that third party tools are not suitable.
Granted, you will need to evaluate that k-mer profile similarity has some relevance to your intuitive notion of "similarity", which it may not. It would clarify somewhat to know more about your desired application.

After more search, I go for blastp-short with MHC peptides as query and my full-length proteins as database. Please share if you think it has problems.


text:
How to extract all variant alleles that do not match "./." from the GT column of a vcf file?
For the two vcf files linked below, I cannot find any variants in the GT column other than "./.". Is it possible to confirm if the GT column of vcf files have been annotated (i.e variants listed as "./1", "1/." or "1/1")?
vcf files: (https://zenodo.org/records/6558593, LUD.TH179.PASS.dbsnp_cosmic.vcf.gz and LUD.TH238.PASS.dbsnp_cosmic.vcf.gz).
here are a few terminal commands I have tried. The output for the commands below are all 0, with no variant information under the columns of the outputted vcf data. I would assume, it means that the GT column was not annotated, but I am not sure.

remove homozygous ref genotype from multi-sample vcf

bcftools view -i 'GT[*]="alt"' LUD.TH179.PASS.dbsnp_cosmic.vcf.gz | less -SN

extracting heterozygous snp from a vcf file

vcftools --gzvcf LUD.TH179.PASS.dbsnp_cosmic.vcf.gz --extract-FORMAT-info GT | grep "0/1"

ID heterozygous variants in VCF file using vcftools

vcftools --gzvcf LUD.TH179.PASS.dbsnp_cosmic.vcf.gz --het
output (cat out.het)

INDV O(HOM) E(HOM) N_SITES F (EMPTY)

Can you clarify the issue please? Yes, it looks like this vcf only has ungenotyped variants. What exactly are you asking?

@coderQ please carefully consider the answer below: I think you'll find it works, then please consider upvoting and marking as "accepted" (it helps everyone).

I checked one of your files, LUD.TH179.PASS.dbsnp_cosmic.vcf.gz, and can confirm that none of the samples listed in the VCF have a genotype for any of the variants. This isn't surprising, it is often hard if not impossible to genotype variants from somatic samples since you cannot know if the, for example, 10% allelic balance you see means that this is homozygous in 10% of the cells or heterozygous in 5% or any other combination of the two.
In any case, this is how I confirmed that no variants were genotyped. First select all lines that are not headers (do not start with #), then fields 10 to the end of the line, convert tabs to newline characters (this leaves me with the sample fields, e.g. ./.:.:.:.:.:.:.), then cut at the first : to get the genotype, sort and count:
$ zgrep -v '^#' LUD.TH179.PASS.dbsnp_cosmic.vcf.gz | cut -f 10- |
   tr '\t' '\n' | cut -d: -f1 | sort | uniq -c
537123250 ./.

As you can see above, the result was 537123250 sample/variant combinations, all of whose genotype was ./..


text:
Multiple genome assemblies of the same bacterial species
I have some RNA-seq data where there are reads from the "host" as well as from several bacteria species. In this experimental context, I am interested in the host associated reads and the latter is regarded as "contaminants".
I am using BBSplit to split host- and bacteria-associated reads using my genome of interest as well as a "repo" containing thousands of bacterial genomes that are likely to be present in the host.
This repo has multiple genome assemblies for the same species making the BBSplit step computationally expensive. The total size of the repo is approximatly 25GB, more than 25X the human genome! I am wondering what would be the best criteria to pick one assembly of a particular genome over another assembly.
All in all, different genome assemblies of a given species have more or less similar total lengths, %GC, number of contigs meaning that picking any would probably be fine but not being a microbiologist, I would like to hear expert's opinions on this. My current take is to take the assembly with the largest size for each species.

Have you tried testing a couple of randomly chosen ones? Do your results vary significantly?

This is rather difficult to test as this would require comparing count matrices generated after removing reads associated with one bacterial genome assembly for another. There are hundreds of species and each has multiple assemblies. The biggest problem is the computational overhead with the mapping with BBSplit. I believe the best way to proceed would be generating a pan-genome for each species and then map to these but I don't have much experince with working with pan-genomes.

do you have the genome of your organism of interest/host? Or even a transcriptome of a related one? You may also cluster reads in your FASTQ files using clumpify from BBMap prior to any mapping/species assignment to speed upthe subsequent steps.

Fair enough. I was just thinking that if you tried the analysis with one genome and then with another and compared the results, it might give you a vague indication of whether it's worth playing with the many genomes or if, for your purposes, any of them will be good enough.

"My current take is to take the assembly with the largest size for each species."  .... Well yes possibly for bacteria - its complicated ..... We'd need more biological context and what the central biological problem is. Basically whats the host? Is this an infection?

Thanks @darked89, I was not aware of clumpify, I will take a look. And my host organism is human.

@haci imho you can just map all the reads to the human genome using splice aware mapper (i.e. STAR). The overwhelming majority of reads from non-eukariots will simply not map to a human genome. If you are not interested in them, there is no need to perform a costly filtering.

Hi @darked89, while your suggestion would be OK for most of the genes, there might and will be genes that will be affected. At least I cannot claim otherwise. As a matter of fact, an initial quick and dirty analysis shows that my count matrices generated with our w/o filtering has a correlation coef of 0.93, much less than I anticipated. I will revive this thread once my colleagues compare the DEG lists to see the "real effect".

@haci what is your read length? Also: what are genes for which there is such cross talk between human and bacteria? Last but not least: unmapped reads do not in influence your human read counts => you may speed up filtering of the putative bacterial reads using just the reads which mapped to human genome. Bonus point: this will let you check if the filtered out reads really are bacterial or human.

@darked89, the whole point of this exercise is that I don't have to "worry about" such a cross-talk, if any. I know that unmapped reads would not have an effect, my problem is the opposite, "the origin of reads that map to the human genome". I am already performing what you suggest and I can tell a simple matrix subtraction of count tables generated w or w/o filtering does not give a matrix full of zeros or sth similar as I was initially expecting, so there is some effect. And I have 50bp SE reads, that is another problem, I would have felt much safer with longer PE reads.

I tend to agree with @darked89 here. Most reads that come from bacteria will simply not map when you map to the host. Moreover, for those reads that actually would map to both, host and bacteria, you have no way of knowing where the respective RNA actually originates from. Removing those would assume that they certainly not from the host. The fact that some reads that map to human also map to a huge collection of bacteria does not really prove that they are of bacterial origin.

A couple of ideas:

Use kraken2 to classify the reads. If you are able to convert your "repo" to kraken2 format, then the matching will be very fast and relatively accurate. Also, using bwa-mem on a joint database could work.
I am not a big fan of custom and not well defined "repos". Why don't you just use the official RefSeq repository of prokaryote genomes? (https://www.ncbi.nlm.nih.gov/refseq/about/prokaryotes/). I don't remember its size, but I remember using it without any problem.

Thank you for the suggestions. 1) I have used Kraken2 as well as Centrifuge for read classification. Apparently the latter is better for shorter reads. Anyway, both failed to call more or less half the reads, probably because my read size of 50bp. 2) What I am using is the human oral microbiome database, they fetch genomes from RefSeq. And then again, I don't want to map to everything but to a collection of "good representatives".

Some more "random" ideas, in addition to the idea of choosing the largest assembly for each species:

you can select for each species the Reference assembly, which is indicated in RefSeq.
I wouldn't be against darked89 suggestion (map to humans only) but I also understand that this should be done only in case of emergency.
Have you considered mapping with STAR on the joint human- bacterial database?

Thanks! Taking the reference for each species is definetly a good idea. And reagerding your third point, can you please elaborate, where can I find this human-bacterial database? Good old Google failed me this time. For STAR, I would need a genome FASTA and associated annotation GTF file and I belive such a combined genome would be GBs big. And regarding your 2nd point, how do you explain the fact that I end up with different count matrices with and w/o filtering for bacterial reads?

Regarding 3rd point: you merge human and all the bacterial fasta in a multifasta and the build the STAR index with that. You are right it will be several gigabytes, but it is possible that once you create the index (IF creating the index gives no problems) STAR will be fast enough in mapping. Regarding 2nd point: I think you are right in supposing that some gene portions may be not too dissimilar from bacterial transcripts, and thus it would be better to map on all the genomes at the same time.

My mapping attempts with BBMap took days with a the full set of genomes (hence the question in the first place) not to mention many problems I had encountered while creating the BBMap index. I don't think doing it with STAR will be much much faster. Thanks anyway.

What I would suggest, without knowing much biology about the system is it's a tree problem. The difficulty is knowing whether the repo is prokaryotic. To solve it I'd take the 16S gene and make a tree - each 16S sequence contains its genome ID. I'd then draw a tree, e.g. IQ-TREE is good and fast. This will describe the overall diversity and provides and informed basis of how to represent each species.

Your approach certainly has merit because again if these are bacteria they go through gene loss - which can cause problems. Selecting the largest genome ensures you minimise this loss.

What I would do is combine 1 and 2. This is to label the taxa that has the largest genome and then produce the 16S tree. It will then be obvious by manual inspection if your approach is representative of the genetic diversity. If it isn't .... you can do the whole thing manually (could be painful), the other approach is to use a tree reading package/module/library. I use ETE3, you'd use ape (@haci's from memory is a big R codes). The issue is traversing the tree from the 'biggest genome' to the most distant genetically related member of the species. What you do is identify the MRCA (most recent common ancestor) for that species and pick a taxa that is in the opposite clade to the biggest genome and any member will do. Thus this is automating the process via Python or R code so thousands of genomes/species can be screened.
To try and be clearer, what you are doing in the code is moving the tree one node at a time and checking if all the members below it are the same species. As soon as the answer is 'no', it's the previous node thats the MRCA for the species. In a bifurcating tree that divides the data into its two biggest groups of genetic diversity.
This approach ensures you get a reasonable representation of the species, by selecting two members of it. One is the largest genome.
It's a while since I've used ETE3, so I can't remember the methods, there's probably a MRCA method for a given species. You can do this empirically, doing a recursive climb up the tree from the biggest genome until your monophyly contains a difference species. The previous traversal is the MRCA, then you ask it to divide the species into two clades and any member of the clade without the biggest clade is what you want.
I agree that it needs a bit of "tree-thinking" to do this ... but it's doable. You then use the genome IDs to filter your repo and produce a greatly slimmed down repo.
You could alternative switch to an HPC method (e.g. AWS) using the entire repo
The alternative approach is know a lot about the biology/ pathogenesis of the problem.

Well the answer is correct, regardless of the number of downvotes. Explaining it clearly - particularly in code - is a lot of text. Basically to automate a representative genetic diversity data set is not trivial, unless there's an external biological rationale, but here nothing specified.

Thanks @M_ for the detailed answer. The contaminants are bacteria and the host is human. I think I will first try your suggestion of using 16 trees to see if taking the largest genome was a good call. The fact that I have over 500 species, I somehow need to automate the representative/not-a-good-representative call.

BTW, regarding the downvote in this answer, I don't really find it useful downvoting but then not explaining why. In this case, here is a complex answer in a domain that I am not experienced (hence the question) and now I don't not know which part might be "off". Giving a reason for the downvote would spark a discussion that the community would benefit.