# Question
Title: Bedtools genomecov gives me coverage files showing zeros at all nucleotide positions
Also posted on biostars

I have bam files that I've obtained through STAR alignment. Here is the first 10 lines of one of my bam files:

```
GWNJ-1013:671:GW2304233104th:2:1126:15501:17879 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1126:15130:18427 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1206:9272:29857  0       chromosome      1       255     28S108M14S      *       0       0       GGCGGTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGAGATCGGAAGAGCG  FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:106        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:30734:10692 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:31412:33411 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFF::FFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFF:FFFF,FFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFF:FFF::FFFFFFFF,FFFFFFFF:FFFFFFFFFFF:F,FFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1330:21938:10708 0       chromosome      1       255     16S134M *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,FFFFFFFFFFFFFFFFF:F  NH:i:1  HI:i:1  AS:i:132        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:17074:4726  0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:11939:11584 0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1372:24740:30843 0       chromosome      1       255     16S48M  *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAA        FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        NH:i:1  HI:i:1  AS:i:47 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1401:11107:19319 0       chromosome      1       255     5S144M1S        *       0       0       ATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGGAATTGA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF  NH:i:1  HI:i:1  AS:i:142        nM:i:0

```

This is the script I'm using that contains the `genomecov -d -pc -ibam` command:

```
BAM_DIR='./sorted_bam_files'

# Loop over each BAM file in the directory
for file in "$BAM_DIR"/*_high_mapq_reads.sorted.bam
do
    echo "Calculating coverage for $file..."
    # Get the basename of the file
    base=$(basename "$file")
    # Calculate coverage
    bedtools genomecov -d -pc -ibam $file > "${base%.sorted.bam}.coverage"
done

```

This is the output I keep getting (first 10 lines only shown). I've tried this with at least 30 other bam files but the output is always the same.

```
chromosome      1       0
chromosome      2       0
chromosome      3       0
chromosome      4       0
chromosome      5       0
chromosome      6       0
chromosome      7       0
chromosome      8       0
chromosome      9       0
chromosome      10      0

```

Any ideas on why this keeps happening? I don't think there's something wrong with my bam files because I was able to use them for DGE analyses and also checked to see if most reads were getting unmapped but they weren't. Any help is appreciated.

# Answer
> 2 votes
option `-pc` is for

```
-pc     Calculate coverage of pair-end fragments.

```

but your reads are single-end reads

---
Tags: rna-seq, bedtools
---


# Question
Title: How to prove correlation between gene expression and functions using omics data in bioinformatics?
I am new to bioinformatics, coming from a background in mathematical biology. I am currently participating in a project where the Principal Investigator (PI), who is a biologist, has tasked me with proving a correlation between the expression and functions of two or more elements. Specifically, one of them is a kinase enzyme and the other is mTOR. The PI has asked me to use omics data or databases for this purpose.

As someone with limited experience in bioinformatics, I am struggling to understand how a biological pathway or correlation can be definitively proven, beyond basic statistical methods. Additionally, the ultimate goal of this research is to develop treatments using substances that our lab has discovered, which have shown promise in killing specific cancer cells in murine and cell lines (although our data on this is limited).

My question goes as followes- is it possible to demonstrate a connection between these elements using only databases and omics methods? Why might a dynamic modeling approach, such as ordinary differential equations (ODEs) not be enough?

I am feeling quite confused and would greatly appreciate any guidance or direction on how to approach this problem effectively.

Thank you.

# Answer
In general terms data mining is needed, then supervised learning and finally ODE, because you are performing DE against enzyme kinetics and thats doable, but you'll need a lot of biochemical knowledge.

I would first take the 'omics' data sets and perform a PCA analysis on it and see where your mTOR and kinase genes are located. PCA is known as unsupervised learning and can be accessed by Seurat or better using Scanpy packages. That will help establish the overall behaviour of your key enzymes pre- and post- treatment ... against the control.

Basically you want:

* Data sets of cancer mice - pre-treatement
* Data set of cancer mice under treatment (drugs targeting kinases)
* Healthy controls

See what unsupervised learning looks like between the 3 groups and where your enzymes are. The PCA plots then compared for the 3 groups and you want to see the targets you are interested in changing groups, especially the treatment group shifting between the control and the diseased unsupervised plots. This would provide a basis for exploring supervised learning (outlined below).

---

**Supervised learning** What you will do is look at gene expression and train it against your phenotype of interest (cancer and recovery from cancer). The training 'target' is the drug the mice are being treated with for each given data set, obviously they'll be a "untreated" and "control" categories. If your model trains correctly \>0.85 to 0.9 accuracy ... post-back and I'll explain how to extract the answer. This approach will extract the biochemical pathway for how the cancer treatment works, but requires more work than described here and that provides the framework for the DE.

Again, a lot of serious knowledge of biochemistry is needed to deploy ODE, unless you've got a brilliant *a priori* \- then you could dive straight in. To do that you'd need to know the exact mechanism of how these enzymes are involved in cancer and then parse that data from your data sets (totally different approach to what is described above). The problem with this approach is if someone knows both the maths and genomics (like a reviewer or examiner) and you've made an error, that will that will be a problem. If you've been through an orthodox approach (above), its a lot stronger justification for the ODE, thus difficult to criticise.

**Note** The key is what data set types comprise 'omics', I'm assuming this is RNA-based (thats normal), but it might not be.
> 2 votes

---
Tags: genome, cancer, pathway, correlation, modelling
---


# Question
Title: Extracting P-Values for Co-Expressed Genes Within WGCNA Modules
I'm facing a challenge regarding the extraction of co-expressing genes from each module. After identifying modules, I've been using exportNetworkToCytoscape to visualize genes that co-express together. However, I solely seek the p-values of the co-expressed genes in each module, similar to the format below:

```
fromNode        toNode      corr-value      p-value
ENSG234        ENSG567       0.90           0.001
ENSG987        ENSG456       0.55           0.03

```

While the exported data provides useful information such as "fromNode," "toNode," "weight", it lacks the p-values (more important) and corr-value I require for my analysis.

I attempted to use tomSimilarityFromExpr with 'corType = "pearson"' to calculate pairwise correlations between genes, but I'm uncertain about the subsequent steps to obtain the desired p-values. My goal is to export the co-expression genes of each module, as facilitated by exportNetworkToCytoscape, while also having the corresponding p-values for these co-expressed genes within each module. I'm not interested in linking the modules to external variables; rather, I solely seek the p-values of the co-expressed genes in each module.

Even if there is another way to obtain what I want without going through tomSimilarityFromExpr and exportNetworkToCytoscape, this suggestion would be welcome!

Could you kindly offer guidance on how to proceed in obtaining this information? Thanks in advance!

# Answer
> 2 votes
The simplest way is to calculate the correlation matrix of genes in the module you want and use that as input to `exportNetworkToCytoscape`. Something like this:

```
inMod = moduleColors=="blue"  ## replace 'blue' with the color or module label you want
modExpr = datExpr[, inMod];
modCor = cor(modExpr)
cytoData = exportNetworkToCytoscape(modCor, threshold = 0, ...)
## Replace ... above with appropriate arguments

```

Then add the p-values to the edge data frame returned within `cytoData` using function `corPvalueStudent`, for example as

```
edgeDF = cytoData$edgeData
nSamples = nrow(datExpr);
edgeDF$pValue = corPvalueStudent(edgeDF$weight, nSamples)

```

You can then manipulate the correlations and p-values further as necessary, e.g. export them into a csv file.

---
Tags: r, wgcna
---


# Question
Title: How to know if FASTQ/BAM is from reference genome (FASTA)?
I'm new to bioinformatics. I have a problem in which I have a FASTA reference genome and lots of reads in FASTQ files. Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome. I also have access to the BAM/BAM.BAI files associated with each FASTQ, although I can't know how they were aligned.

Thanks

# Answer
> 3 votes
Your post has multiple questions but one of them is not really a question so I'll address everything here:

1. "Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome"

Just align the reads. Contaminants should not align well and should be discarded owing to poor mapping quality score.

2. How to figure out if the FASTA file was used for the alignment

This is not really a 100% perfect solution but the BAM header should contain the exact command used for the alignment, which should contain the ref FASTA with path. You can also compare contig names and lengths between the FASTA and the BAM `@SQ` lines.

# Answer
> 2 votes
In general, this answer wouldn't be more correct than @Ram RS' answer but here is how I would approach:

* I would start everything from scratch, at the end of the day you will most probably need to document the exact genome version and its associated annotation file. It looks like you do not have this information at present.
* I don't know the experimental set-up but if you are worried about "contaminants", you can have an extra alignment step using `BBSplit` (link), it is capable of mapping to multiple genomes simultaneously.
* I would give a try with `nfcore/rnaseq`, it abstracts away a lot of complex bioinformatics steps so that you can focus on understanding what these tools do instead of dealing with installing and running these tools. For example you can just use the `--bbsplit-fasta-list` argument to get rid of your contaminants.

---
Tags: sequence-alignment, fasta, fastq, read-mapping
---


title: Gene list from eclip-seq .bed files

text:
# Question
Title: Gene list from eclip-seq .bed files
I have some eclipseq data, three IP and three input datasets, and have the alignments for all three. I used Bedtools "Intersect Interval" with the genome gtf file (hg19) I used to align them with, and that creates a bed file with gene names for each peak. I then wrote an R script to sum the clipper score for each gene within a sample, then do a simple student's t test across groups. I now have a gene list and p values for differential expression, but not sure if this is correct.

I've tried looking at a lot of sources but cant really find any canonical way to generate a gene list from bed files.

I've also tried IDR and it doesn't work on the clipper output, unfortunately.

Thanks!

Edit: eCLIP-sequencing, it's essentially a variant of iCLIP-sequencing and data analysis methods are transferrable, to my knowledge.

I have alignments (bam files) made with Star aligner against Hg19. After deduplication and sorting, I called peaks using Clipper (a peak-caller made by the creators of eCLIP). It outputs bed files with scores for each sample.

My question is how exactly do I turn this bed file with scores into a gene list? Is my method above correct?

# Answer
> 1 votes
As someone who works with eCLIP data from ENCODE, I think I can help out.

You need to think of this more similar to ATAC or RNA-Seq, whereby you actually get a counts matrix, normalize that, and then apply math.

Another good idea would be to look through this eCLIP analysis document provided by Gene Yeo of ENCODE for more specific ideas.

https://www.encodeproject.org/documents/1f171ac6-a36a-41ac-b632-741aeb47aad2/@@download/attachment/eCLIP\_analysisSOP\_v2.3.pdf

# Answer
> 0 votes
As others have said it's not very clear what you're asking here, but I'm assuming you want to know how to determine relevant genes that were enriched in your IP.

As per the group behind eCLIP, Van Nostrand et al. 2020 : "standard peak analysis used the set of peaks identified as irreproducible discovery rate (IDR) reproducible and meeting fold-enrichment (≥ 8-fold) and significance (p value ≤10−3) in immunoprecipitation versus paired size-matched input."

If those cutoffs are too stringent you could change them to fit the context of your project. And what do you mean by differential expression? It sounds like you're addressing binding.

---
Tags: peak-calling, fold-change
---


# Question
Title: 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.

# Answer
> 2 votes
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:

---
Tags: r, upsetr
---


# Question
Title: 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?

# Answer
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.
> 3 votes

---
Tags: sequence-analysis, deseq2, gene-expression, go-enrichment
---


# Question
Title: Finding annotated counterpart after BLASTn with efetch (Biopython)
I am creating a pipeline for the identification of unknown transcripts. After a local BLASTn search of the transcripts, I have a large list of the respective hits with different genomes. I have the accession code of the respective genome, start/end coordinates, and the strand where the hit is found. Now, I would like to find if the genome where the hit was found is annotated within the coordinates returned by BLAST.

I have found the possibility of using efetch to search for each genome via Biopython, but this only returns the whole annotation, even though I limit the search to the specific coordinates with:

```
Entrez.efetch(db="nuccore", id="CP054847.1", strand="2", seq_start="326053", seq_end="326786", retmode="xml")

```

Of course, I could parse the XML and look for possible annotations within the given range, but not only would this be quite computationally intense since I expect around 10.000+ hits to be searched, but it also seems to return the same results for the plus/minus strand, even though the annotations should be strand-specific.

I have also intended to use the standalone tool of efetch, but it would be easier to integrate the search to the pipeline using Biopython.

Do you have any ideas why `Entrez.efetch` is returning the whole `xml` or is there any fast way of filtering the output? Do you have maybe another alternative/idea for this step? I would appreciate any suggestion!

# Answer
> 1 votes
It seems like the `seq_start` and `seq_end` are sensitive to the direction. Try swapping the two.

# Answer
> 1 votes
Lets assume you can get the `blastn` output in a `BED` format:

```
genomeA_contig001    hit_start    hit_end    query_transcript    score   strand
genomeB_contig021    hit_start    hit_end    query_transcript    score   strand

```

which will contain all the unknown hits to all genomes.

All you need is a set of `GTF` files specific for `genomeA`, `genomeB` and so on where you add the suitable prefix (`genomeA_` for `GTF` describing genes/exons in that genome and so on). For transcript hits you may start with just the annotated exons. Next run:

`bedtools intersect -wao -a all.transctipts.hits.bed -b all_annotated_exons.gtf my_result.tsv`

to get all the overlaps with number of overlapping bases. Filter and sort the `my_result.tsv` possibly importing it into a data frame using `pandas` or `polars`.

---
Tags: python, sequence-annotation, blastn, entrez, efetch
---


title: How to parse hmmsearch output?

text:
# Question
Title: 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

# Answer
> 2 votes
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

# Answer
> 1 votes
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

```

---
Tags: linux, biopython, functional-annotation, parsing, hmmer
---


# Question
Title: 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)
```
</code></pre>

# Answer
> 3 votes
```
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.

---
Tags: rna-alignment, star
---


# Question
Title: Advise on building an effect ML model for predicting important proteins for drug response
I want to create a model to predict proteins which could be associated with drug response in cancer cell lines. I have cell line proteomics data, with compound screening data they have gained for a number of these cells. My aim is to relate these proteins to drug response. I am using this exercise as an introduction in the use of ML in biology.

To this end, I have created dictionaries, with selected compounds as the keys, and their respective datasets. Cells are rows, proteins are columns, and the drug response is the `y` variable.

I have approached the task as a regression task. The aim being to predict the ln(IC50) - a measure of the concentration required to induce cell death - and subsequently pull out the features (proteins) which have made the largest contribution to this prediction. As my baseline model, I have attempted linear regression. However, my model performance has consistently been very poor. A problem with this dataset is the large number of missing values (heatmap below - yellow = missing). However, even when trying a simple feature selection strategy with sklearn (first function), the explained variance is around 0.3 for the best performing compound. I have plotted the residuals in attempt to dig into this further, but the spread looks quite distributed.

Any advise/suggestions would be appreciated. Or any further reading/papers relating to ML in biology which could help a beginner build and use effective models.

Training and test functions shown below.

```
    def process_training_and_test_data_for_LG(df):
        training_datasets = {}
        test_datasets = {}
        for drug in df.keys():
            df[drug].dropna(axis=1,inplace=True)
            print("at "+drug)
            X = df[drug].drop(["cell_line_name", "ln_IC50","putative_target"], axis=1)
            y = df[drug]["ln_IC50"]
            X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

            feature_selector = SelectKBest(score_func=f_regression, k=50)
            feature_selector.fit(X_train, y_train)

            X_train = feature_selector.transform(X_train)
            X_test = feature_selector.transform(X_test)

            #Add data to relevant dictionaries
            training_datasets[drug] = [X_train, y_train]
            test_datasets[drug] = [X_test, y_test]
        return training_datasets, test_datasets

def train_models(training_data):
    trained_models = {}
    for drug in training_data.keys():
        #initiate model
        lm = LinearRegression()
        #extract drug specific training data
        X_train = training_data[drug][0]
        y_train = training_data[drug][1]
        #fit model
        model = lm.fit(X_train, y_train)
        #Append dict of trained models
        trained_models[drug] = model
        print("Model for {} trained".format(drug))
    return trained_models

def test_models(model_dict, test_data):
    prediction_dict = {}
    for drug in model_dict.keys():
        drug_model = model_dict[drug]
        X_test = test_data[drug][0]
        predictions = drug_model.predict(X_test)
        prediction_dict[drug] = predictions
    return prediction_dict

```

# Answer
> 1 votes
The missing values might be an issue indeed. You might want to use imputation methods, e.g.

```
from sklearn.impute import KNNImputer
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer

def impute_df(df, feature_columns_to_apply, strategy):
    if strategy == 'mean':
        for col in feature_columns_to_apply:
            df.loc[:, col] = df.loc[:, col].fillna(df[col].mean())
    elif strategy == 'median':
        for col in feature_columns_to_apply:
            df.loc[:, col] = df.loc[:, col].fillna(df[col].median())
    elif strategy == 'knn':
        knn_imputer = KNNImputer(n_neighbors=5)
        df_tmp = pd.DataFrame(knn_imputer.fit_transform(df[feature_columns_to_apply]), columns=feature_columns_to_apply)
        df.loc[:, feature_columns_to_apply] = df_tmp[feature_columns_to_apply]
    elif strategy == 'iterative':
        iterative_imputer = IterativeImputer()
        df_tmp = pd.DataFrame(iterative_imputer.fit_transform(df[feature_columns_to_apply]), columns=feature_columns_to_apply)
        df.loc[:, feature_columns_to_apply] = df_tmp[feature_columns_to_apply]
    return df

df_all = impute_df(df_all, feature_columns_to_apply, 'mean')

```

(IterativeImputer is heavy and slow) These are just sample strategies, you might want to use other methods, see e.g. here.
Also consider other models such as random forest / XGBoost, they outperform linear regression in many cases.

# Answer
> 0 votes
The gist looks fine ..

* `pandas` \- looks fine
* `SelectKBest` \- cool
* transformation (obviously there's a lot of them)
* standard ML

A common ML situation is, the classifier working but the accuracy is notably lower than expected. My personal *first thing to check* is the use of *linear regression* as a baseline model. The scatterplots give a gist, we don't know in particular what data points belong to what classification - thats crucial. What's helpful are coloured dots to assess how they sit against the regression line.

What I think has happened is there's been a PCA via `SelectKBest`, unsupervised learning, which doesn't account for much of variance (*cf* 0.3), i.e. this is not an accuracy score.

I've three suggestions:

1/ An acceptable baseline model often used is `k-nearest neighbour` ..

```
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
clf = KNeighborsClassifier(n_neighbors=3)
clf.fit(X_train, y_train)
print("Test set accuracy: {:.2f}".format(clf.score(X_test, y_test)))

```

I think this will produce much better results with an accuracy of \>0.7

2/ However, I would check the following (please see **caveat**),

```
for drug in df.keys():
    df ....
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

```

Just make sure that X\_train, etc ... is supposed within the loop, rather than immediately outside it (four spaces to the left). Looking at your outputs they look fine.

The thing to check is its not:

```
for drug in df.keys():
    df....
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.33, random_state=42)

```

The above loop makes sense to me.

The original loop position `train_test_split` looks like the `key` loop might repeatedly over-ride the `X_train, X_test, y_train, y_test` outputs.

3/ Very final point `SelectKBest` followed by `fit`. What about `fit_transform` instead? If the issue is the features are resulting in low description of the variation, maybe thats key here?

---

**Caveat** I do stress that its hard to comment beyond general terms because a good biological understanding is needed of the goal of the ML training. In context its not my area and I haven't read any of the papers.

---
Tags: python, machine-learning, proteomics, modelling, linear-regression
---


# Question
Title: 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.

# Answer
> 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)

```
> 0 votes

---
Tags: python, machine-learning, algorithms
---


title: How to debug a bash script

text:
# Question
Title: 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

# Answer
> 4 votes
> 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

```

---
Tags: vcf, bash, linux, wgs, for-loop
---


# Question
Title: 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.

1. remove homozygous ref genotype from multi-sample vcf

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

2. extracting heterozygous snp from a vcf file

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

3. 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)`

# Answer
> 2 votes
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 `./.`.

---
Tags: vcf, genomics, bcftools, vcftools
---


# Question
Title: Making a volcano plot look less cluttered
I have generated a volcano plot with a differential expression file.

Code for inputing file:

```
library(tidyr)
final_tumor <- df %>% pivot_longer(., -c(Feature.ID, Feature.Name), names_to = c("set",".value"), names_pattern = "(.+)_(.+)")

# A tibble: 80 x 6
Feature.ID Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>      <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 a          A            Cluster.1    0.000961            0.292           1
2 a          A            Cluster.2    0.000902            0.793           1
3 a          A            Cluster.3    0.00181             1.46            0.758
4 a          A            Cluster.4    0.000642            0.269           1
5 b          B            Cluster.1    0.000320            1.95            0.910
6 b          B            Cluster.2    0.00180             4.77            0.154
7 b          B            Cluster.3    0                   2.19            1
8 b          B            Cluster.4    0                   1.66            1
9 c          C            Cluster.1    0.00128            -2.01            0.0467
10 c          C            Cluster.2    0.00632             0.352           1
# … with 70 more rows

```

Here the output of head(final\_tumor)

```
# A tibble: 6 x 6
Feature.ID        Feature.Name set       Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>           <dbl>            <dbl>            <dbl>
  1 ENSG00000227232.5 WASH7P       Cluster.1     0                   1.50            1
2 ENSG00000227232.5 WASH7P       Cluster.2     0                   1.73            1
3 ENSG00000227232.5 WASH7P       Cluster.3     0                   1.77            1
4 ENSG00000227232.5 WASH7P       Cluster.4     0.00114             4.30            0.293
5 ENSG00000227232.5 WASH7P       Cluster.5     0                   2.15            1
6 ENSG00000227232.5 WASH7P       Cluster.6     0                   1.22            1

```

And here is the output of tail(final\_tumor):

```
# A tibble: 6 x 6
Feature.ID        Feature.Name set        Mean.Counts Log2.fold.change Adjusted.p.value
<fct>             <fct>        <chr>            <dbl>            <dbl>            <dbl>
  1 ENSG00000210196.2 MT-TP        Cluster.6       0.0699          -0.202           0.790
2 ENSG00000210196.2 MT-TP        Cluster.7       0.0801           0.0386          1
3 ENSG00000210196.2 MT-TP        Cluster.8       0.0711           0.0875          1
4 ENSG00000210196.2 MT-TP        Cluster.9       0.0152          -2.31            0.00127
5 ENSG00000210196.2 MT-TP        Cluster.10      0.0147          -2.30            0.00612
6 ENSG00000210196.2 MT-TP        Cluster.11      0.122            0.762           1

```

Here is the code for generating volcano plot:

```
library(ggplot2)
library(ggrepel)
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.05),
                  aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"),
                  aes(label = Feature.Name), color = "red")

```

Generating the ggplot with all the data was taking forever, it was suggested to select random rows of the dataframe:

```
rep_final = final_tumor[sample(1:nrow(final_tumor), size = 1000), ]

```

Then I did:

```
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.05), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

```

Everything is very tight and cluttered. I want to show the top differentially expressed genes. How do I go about simplifying the plot?

Is there a way to do this on EnhancedVolcano, pulling out a certain gene?

```
EnhancedVolcano(final_tumor, lab = as.character(final_tumor$FeatureName), x = 'Log2.fold.change', y = 'Adjusted.p.value', xlim = c(-8,8), title = 'Tumor', pCutoff = 10e-5, FCcutoff = 1.5, pointSize = 3.0, labSize = 3.0)

```

# Answer
**The element that is clustering your plot the most is the ggrepel labels**, so that's what you need to work on.

To get to a readable number of labels, **I would not remove random genes**: you risk losing important info, be it genes that could be of interest to you or some of the most regulated. Incidentally that seems to happen on your example: you try plotting in red the label for Casp14, but it does not show on the plot, probably because the gene was removed in your sampling.

What you can do however is, as dc37 suggests, **increase the thresholds for the labels** (you can probably play around with the thresholds to see when it does become readable). This should also help with the plotting: you mention that generating the plot was taking a long time, but with the number of genes you are looking at (~10k-30k?), the geom\_points() should not be a real problem (maybe a few seconds), however the ggrepel is the most computationnally expensive so that's what you should try to reduce.

This should get you started:

```
ggplot(final_tumor, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(final_tumor, Adjusted.p.value < 0.001 & abs(Log2.fold.change) > 1 ), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name), color = "red")

```

Notice everything plots the "final\_tumor" object and not the "rep\_final". Feel free to be even more stringent with the cutoffs if it's still taking too long to plot.

A bit outside of the score of your question, but I also wanted to mention that you should never be trying to label too many genes on a volcano plot: this type of graph is useful to visualize your data's p-value/FC distribution and highlight a subset of genes in this context (most regulated, genes regulated in another experiment, pathway of interest...), but if you label too many of them it will never be readable.

If you really want to be able to indentify all the genes on your graph, I would try an interactive visualization (like with plotly). I'm not going to detail the code here but you could look in that direction! (However, as with ggrepel, you might need to remove some genes so that it does not "take forever").

Have fun :)
> 2 votes

# Answer
Aside from the `ggrepel` labels, another reason your volcano plot looks cluttered is the presence of the 4 very highly significant genes on the top right, which lead to the rest of the genes being squeezed into the x axis. For visualization purposes, you could try applying an inverse hyperbolic sine transformation (`asinh()` in R) to both your x/y axes, which should alleviate this squeezing while not distorting too much the overall volcano shape. Just be sure to also apply this transformation to any cutoff lines you are including in the plot.
> 2 votes

# Answer
As we discussed in comments in your previous post: Pulling out a certain gene in a volcano plot, in order to have less labeled genes on your volcano plot, you can specify more stringent cutoff.

For example, you can decide to display only genes with an absolute fold change superior to 2 by selecting the appropriate subset of your data:

```
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

```

Alternatively, you can maybe decide that a p value inferior to 0.05 is not enough stringent and decide to labeled only genes with an adjusted p value inferior to 0.00001 by doing:

```
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, Adjusted.p.value < 0.00001), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

```

Or, you can want to do a combination of both criteria by doing:

```
ggplot(rep_final, aes(x = Log2.fold.change,y = -log10(Adjusted.p.value), label = Feature.Name))+
  geom_point()+
  geom_text_repel(data = subset(rep_final, abs(Log2.fold.change) > 2 & Adjusted.p.value < 0.00001), aes(label = Feature.Name))+
  geom_text_repel(data = subset(final_tumor, Feature.Name == "Casp14"), aes(label = Feature.Name),
                  color = "red", size = 5)

```

At the end, I think it depends of what you are trying to demonstrate by presenting this volcano plot. Do you want to highlight some particular genes? Or do you want to display genes that are the most differentially expressed? If so, how many?

All of this depends of the scientific question underlying your experiments and why you came to make this volcano plot.

Hope this couple of example helps you to figure it out how to customize your volcano plot.
> 1 votes

---
Tags: r, bioconductor, ggplot2
---


# Question
Title: 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.

# Answer
> 3 votes
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.

# Answer
> 2 votes
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.

# Answer
> 1 votes
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.

---
Tags: blastp, peptide
---


# Question
Title: SARS-CoV-2 sequence used in the AlphaSeq Antibody Datasets to predict binding affinity
Currently we are building a sequence based deep-learning model to predict binding affinity between antibody and antigen. For this we are training a sequence based model with **AlphaSeq Antibody dataset** (https://github.com/mit-ll/AlphaSeq\_Antibody\_Dataset). But in this dataset they haven't provided the sequence for SARS-CoV-2 virus.

So I read the paper and there they have mentioned they have tested against **peptide in the HR2 region of the SARS-CoV-2 spike protein** and they have given the sequence as "PDVDLGDISGINAS". (https://www.nature.com/articles/s41597-022-01779-4#Sec3)

But this sequence is very small. I'd like a larger protein sequence that includes non-epitope portions so that the prediction is more meaningful (e.g. predicts binding to the actual epitope when provided a sequence which includes sequence beyond the epitope). I searched online and only found a covid genome sequence which isn't very helpful.

Searching the RCSB PDB website for Spike protein HR2 gave two SARS-CoV-2 results (*RCSB PDB website*).

One of them had sequence like the following:

```
>8CZI_1|Chains A, B, C|Scaffolded Spike protein S2' HR1|Nostoc punctiforme (strain ATCC 29133 / PCC 73102) (63737)
MSHHHHHHGSQTLLRNFGNVYDNPVLLDRSVTAPVTEGFNVVLASFQALYLQYQKHHFVVEGSEFYSLHEFFNESYNQVQDHIHEIGERLDGLGGVPVATFSKLAELTCFEQESEGVYSSRQMVENDLAAEQAIIGVIRRQAAQAESLGDRGTRYLYEKILLKTEERAYHLSHFLAKDSLTLGFAYENQKLIANQFNSAIGKIQDSLSSTASALGKLQDVVNQNAQALNTLVKQLSSNFGAISSVLNDILSRLDKVE

>8CZI_2|Chains D, E, F|Spike protein S2' HR2|Severe acute respiratory syndrome coronavirus 2 (2697049)
KNHTSPDVDLGDISGINASVVNIQKEIDRLNEVAKNLNESLIDLQ

```

So I am not sure which one is the correct sequence for our task. Can someone provide some insights regarding this matter?

Also why are there so many sequences popping up when I search for SARS-CoV-2? Is there a version that can be considered as the root from which other versions are derived?

# Answer
> 1 votes
PDB has a whole page talking about various slices of the covid spike protein. I'd strongly recommend reading that for a better understanding of why people choose different pieces to represent. In that page they include several different solved structures of the protein that may be helpful for your benchmarking.

In terms of the "canonical" sequence, you could look at a reference genome sequence. Here for example is the EMBL reference genome for covid in GenBank format, which includes the protein sequences. I believe in that the spike protein has ID `QHD43416`.

You can additionally look at NCBI, which has a slightly different version of the reference (though they should be functionally identical, with just slightly different identifiers and annotations).

Assuming that you want the ancestral initial isolate of covid, you can find a similar GenBank file here from NCBI.

---
Tags: sequence-analysis, database, covid-19, sars-cov-2
---


title: How many ML tree should I build?

text:
# Question
Title: How many ML tree should I build?
I created a phylogenetic tree of 328 species using RAxML. I run the ML algorithm 20 times and I chose the best tree (the one with the highest likelihood).

The choice of 20 iterations was based on a tutorial that was carrying out the same analysis but on fewer species.

Now I have a doubt that my search is actually not optimal, and could have got better by increasing the number of iterations.

Do you know if there is a way to estimate the optimal number of iterations to find the best robust ML Tree?

# Answer
> 2 votes
Simply use the iterations option (e.g. `-n 1000`) of IQTREE. It's a lot quicker (several logs) and is strongly of identical performance to RAxML.

```
iqtree -s alignment.phy -n 1000 -m GTR+I+G -nt 4

```

It takes phylip, fasta, nexus ... and no-one questions 1000 iterations, so you're safe.

I personally I'd then go on to bootstrap the tree (IQTREE is extremely fast). If you were bootstrapping with iterations RAxML is unworkable, you'd need e.g. 10 iterations per bootstrap (1000 bootstraps) and would be time consuming. IQTREE could easily perform 10 iterations per bootstrap and would accommodate 100 iterations (1000 would be a bit excessive for bootstrapping).

RAxML is an old program now, IQTREE is its successor and massively refactored for maximum likelihood phylogeny.

---

If you don't understand the command line, drop a note and I'll explain the model it's using.

---
Tags: phylogenetics, phylogeny, taxonomy, dna
---


# Question
Title: Proper formatting of .ped files, in order to run ADMIXTURE
I'm trying to run ADMIXTURE on some files I've generated using the Stacks pipeline. ADMIXTURE accepts .bed and .ped files, but when I try to pass the .ped file generated by Stacks to ADMIXTURE (v1.3.0), I get a "PLINK Input file error". I've pasted the head of my (tab separated) .ped file below.

> head populations.plink.ped -n 1
>
> garden QUAC\_G\_SH\_Q1485 0 0 0 0 0 0 T T C C G G T T C C C C G G T

Alternatively, I could convert my .ped file to .bed using plink, and try to pass this to ADMIXTURE--but Stacks didn't provide me with a map file to do this, so I don't think that will work.

If anyone knows what needs to change with my .ped file (remove sample name? population designation?), I'd appreciate any tips. This is my first time running most of these programs.

Thanks!

# Answer
> 1 votes
Admixture requires a '12'-coded .ped file, rather than the allele-coding you're using. 'plink --file \[file prefix\] --recode 12' produced a .ped file that Admixture could use for me.

---
Tags: plink
---


# Question
Title: Problem with Calculating Polygenic Risk Score (PRS) in Plink
I've been struggling with constructing a PRS score. In the pics you can see how my summary stats and my target data look like.

Target Data Summary Stats

I've tried the following code which should be pretty straightforward:

```
plink \
  --bfile target \
  --score summary_stats.txt 1 4 6 header \
  --out final

```

When I run it, I get this message:

```
10499994 variants and 1198 people pass filters and QC. Note: No phenotypes present.

```

The data is binary and the fam file really doesn't have phenotypes (it's -9). Is there any way to still run a PRS?

I'm unsure what I'm doing wrong and would really appreciate some insight! Also, the data has been harmonized already so it should technically work!

# Answer
> 0 votes
I'm a little confused about what you mean by there being no phenotype. To build the score you need a phenotype, although to apply it you might not. (Usually PRS is build with a base dataset and then applied using a target dataset). Here's a good tutorial about doing PRS in plink that might help you:

https://choishingwan.github.io/PRS-Tutorial/plink/

Also this paper is a really good paper that helped me better understand how building PRS works:

https://www.nature.com/articles/s41596-020-0353-1

Hope this helps!

---
Tags: plink, polygenic-risk-score
---


title: arm64 installation of AMRFinderPlus

text:
# Question
Title: arm64 installation of AMRFinderPlus
I was trying to install NCBI AMRFinderPLUS on my Nvidia Jetson using standard Bioconda based installtion. However, all the time it is giving error that the package is not available and I thought I am making some typo or channel is not proper.

```
(amrfinder) jay@jetmimd-desktop:~$ conda install -y -c conda-forge -c bioconda ncbi-amrfinderplus
Channels:
 - conda-forge
 - bioconda
 - defaults
Platform: linux-aarch64
Collecting package metadata (repodata.json): done
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - ncbi-amrfinderplus

Current channels:

  - https://conda.anaconda.org/conda-forge
  - https://conda.anaconda.org/bioconda
  - defaults

To search for alternate channels that may provide the conda package you're
looking for, navigate to

    https://anaconda.org

and use the search bar at the top of the page.

```

Then I remembered about the different architecture. So, now here I am. Is there any way to install this on arm64 architecture?

Thank you for the help in advance!!!!

# Answer
Try the following

```
conda install bioconda::ncbi-amrfinderplus

```

The above does work I tested it (below). There could be an issue with you version of conda (24.1.2 OSX).

Also has the package has been ported for GPU(?) ... *I assume it is, otherwise it's obviously not going to work and conda will not install it*

However, if it ported for GPU and still fails use mamba (I don't use mamba BTW), because it is very slow and that could be where the error is (unusually slow).

My connection initially was unsuccessful but resolved via

> Retrying with flexible solve.

---

To demonstrate I got it working here's the reply ...

```
The following NEW packages will be INSTALLED:

  blast              bioconda/osx-64::blast-2.12.0-pl5262h78c34c6_0
  c-ares             pkgs/main/osx-64::c-ares-1.19.1-h6c40b1e_0
  curl               pkgs/main/osx-64::curl-7.88.1-hdb2fb19_2
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  hmmer              bioconda/osx-64::hmmer-3.3.2-hb19a6fa_3
  libcurl            pkgs/main/osx-64::libcurl-7.88.1-ha585b31_2
  libev              pkgs/main/osx-64::libev-4.33-h9ed2024_1
  libidn2            pkgs/main/osx-64::libidn2-2.3.4-h6c40b1e_0
  libnghttp2         pkgs/main/osx-64::libnghttp2-1.52.0-h1c88b7d_1
  libssh2            pkgs/main/osx-64::libssh2-1.10.0-hdb2fb19_2
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  ncbi-amrfinderplus bioconda/osx-64::ncbi-amrfinderplus-3.11.14-hfee4e81_0
  perl               pkgs/main/osx-64::perl-5.26.2-h4e221da_0
  perl-archive-tar   bioconda/osx-64::perl-archive-tar-2.32-pl526_0
  perl-carp          bioconda/osx-64::perl-carp-1.38-pl526_3
  perl-common-sense  bioconda/osx-64::perl-common-sense-3.74-pl526_2
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.087-pl526h6de7cb9_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.087-pl526h770b8ee_0
  perl-exporter      bioconda/osx-64::perl-exporter-5.72-pl526_1
  perl-exporter-tiny bioconda/osx-64::perl-exporter-tiny-1.002001-pl526_0
  perl-extutils-mak~ bioconda/osx-64::perl-extutils-makemaker-7.36-pl526_1
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.087-pl526h6de7cb9_0
  perl-io-zlib       bioconda/osx-64::perl-io-zlib-1.10-pl526_2
  perl-json          bioconda/osx-64::perl-json-4.02-pl526_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl526h04f5b5a_3
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-0.428-pl526_1
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.428-pl526_0
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl526h1de35cc_1
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.52-pl526h01d97ff_0
  perl-types-serial~ bioconda/osx-64::perl-types-serialiser-1.0-pl526_2
  perl-xsloader      bioconda/osx-64::perl-xsloader-0.24-pl526_0
  wget               pkgs/main/osx-64::wget-1.21.4-ha585b31_1

```

There is a blast gpu and I'd assume the other library/packages have gpu versions. There's a lot of Perl code there BTW, so it's an old package, before NVIDIA was big perhaps?

---

If there is a GPU compliant channel on conda (nanoporetch) - that is impressive: especially given the amount of Perl code. Porting for GPU from CPU isn't trivial.

> Few versions behind

The key thing is whether the upgrades are fixing bugs or refactoring/adding functionality. If it's the latter then I would go for it because the GPU would compensate for the slower runtime (obviously given you don't need the extra functionality).
> 1 votes

---
Tags: software-installation, ncbi
---


title: NEO analysis tutorial and data

text:
# Question
Title: NEO analysis tutorial and data
I was reading Jason Aten's paper Using genetic markers to orient the edges in quantitative trait networks: The NEO software, which states

> The NEO software can be used to orient the edges of gene co-expression networks or quantitative trait networks if the edges can be anchored to genetic marker data. R software tutorials, data, and supplementary material can be downloaded from: http://www.genetics.ucla.edu/labs/horvath/aten/NEO

But it seems the site is down, does anyone know where can I access these data?

# Answer
Like mentioned in the comment above, I would suggest emailing the corresponding author. Scientific articles usually list the email of the corresponding author and this would be the most logical first step in addressing the missing materials.
> 0 votes

---
Tags: wgcna
---


# Question
Title: How to download multiple UniRef clusters with all their members?
We are trying to download multiple UniRef clusters in one batch on the website. However, it seems the members are not downloaded... Could you guide me a way to the solution? Many thanks!

For example, see this link. We would like to download all these clusters including their members. But when selected these clusters, only the representative sequences (one sequence per cluster) were downloaded.

# Answer
The UniProt page apparently only lets you download fasta sequences for the representative protein in each UniRef group:

So you can't get all the fasta sequences from there. However, once you get the list of UniProt IDs, you can get the sequences from those. To get the IDs, first select TSV output:

Then remove everything and add "cluster members" from the "Miscellaneous" section:

This will give you a file that looks like this:

```
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz
Cluster ID  Cluster members
UniRef50_A0A0F6EQC6 A0A0F6EQC6; E3E6D7; A0A3G8RAV1; A0A2K1EFR7; A0A2S6NYI6; A0A7H0Y5M1; A0A2N9YZ66; A0A853Q3J3; A0A8I2C241; A0A0K2FC93
UniRef50_A0A0T6A5I5 A0A0T6A5I5; A0A2V7BZB0; A0A2V6XWF2
UniRef50_A0A2H4TXU1 A0A2H4TXU1
UniRef50_A0A6C2YI31 A0A6C2YI31
UniRef50_B1KXG0 B1KXG0; UPI001F3F1379

```

It's a tab separated file, and the second field is a `;`-separated list of cluster members. Get those into a new file, one per line:

```
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz |
    tail -n +2 | cut -f2 | sed 's/; /\n/g'
A0A0F6EQC6
E3E6D7
A0A3G8RAV1
A0A2K1EFR7
A0A2S6NYI6
A0A7H0Y5M1
A0A2N9YZ66
A0A853Q3J3
A0A8I2C241
A0A0K2FC93
A0A0T6A5I5
A0A2V7BZB0
A0A2V6XWF2
A0A2H4TXU1
A0A6C2YI31
B1KXG0
UPI001F3F1379

```

Save to a file with

```
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz |
    tail -n +2 | cut -f2 | sed 's/; /\n/g > newfile

```

It may be possible for the same protein to belong to multiple clusters, I don't know UniRef well enough to be sure if this is the case. So, to be on the safe side, make the IDs unique:

```
zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz | tail -n +2 | cut -f2 | sed 's/; /\n/g' | sort -u

```

That assumes GNU `sort`, the default on Linux. If yours doesn't have a `-u` option, use this instead:

```
zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz | tail -n +2 | cut -f2 | sed 's/; /\n/g' | sort | uniq

```

Now, you can use these IDs and UniProt's ID mapping tool to get their sequences. Go to the mapping tool and paste in your IDs (this assumes you have no more than 100000 IDs. If you have more, do this in multiple steps.

Click on "Map N IDs" and that will take you to a page like this, where you can now click on "Download" to retrieve the sequences:

The result of all this is a fasta file like this:

```
$ zcat foo.fa.gz | head
>tr|A0A0F6EQC6|A0A0F6EQC6_PAEPO Membrane protease regulatory membrane protein OS=Paenibacillus polymyxa OX=1406 GN=yuaF PE=4 SV=1
METLFWWLLAGGAIFTVVSVLIGDVLGSWLDGALELPGLDWFRPVVLLGAMTAFGGAGVL
LTKYTGLTMGWVVVLALGIALVIGVLVFFAFIKPMANSEVSSGFSMRELTGRIGEVTIPV
PEVGYGEVMIRLGAGNTIHTASSFDHKPLAAGTRIVVVEVAEGVVRVAYLDT
>tr|A0A0K2FC93|A0A0K2FC93_9BACL Membrane protease regulatory membrane protein OS=Paenibacillus peoriae OX=59893 GN=ABE82_18870 PE=4 SV=1
METIFWWLLAGGAIFTVVSVLIGDVLGGWLDGALELPGLDWFRPVVLLGAMTAFGGAGVL
LTKYTGLSMGWVVVLALGIAFVIGILVFFAFIKPMANSEVSSGFSMRELTGRIGEVTIPV
PEVGYGEVMIRLGAGNTIHTASSFDHKPLAAGTRIVVVEVAEGVVRVAYLDT
>tr|A0A0T6A5I5|A0A0T6A5I5_9BACT Membrane protein, putative membrane protein OS=candidate division NC10 bacterium CSP1-5 OX=1640516 GN=XU15_C0006G0101 PE=4 SV=1
MDSFVQLVVGTVVLRPYVVAFLLLYLFAAVADMGWRRTLLFTVVAWSVAFTAESLSIRFG

```

However, I would urge you to think if this is really the best way for what you want to do. If you are trying to find orthologous groups, there are better ways.
> 0 votes

---
Tags: data-download
---