text:
Using CoRAL to reconstruct ecDNA
We are trying to use CoRAL to reconstruct ecDNA. However, we have a problem in CoRAL.py reconstruct, when an error is obtained in infer_breakpoint_graph step. The only difference is reads' article are aligned with winnowmap, meanwhile ours are aligned with pbmm2, but column 3 of both bam files (paper and study) have identical names ("chr5"). However, if reads are aligned with winnowap there is no error. So, I would like to ask if anyone knows what we should modify in order to use aligned reads with pbmm2.
Error message:
Traceback (most recent call last):
  File "/ngs/coral_soft/CoRAL/src/CoRAL.py", line 150, in <module>
    reconstruct_mode(args)
  File "/ngs/coral_soft/CoRAL/src/CoRAL.py", line 23, in reconstruct_mode
    infer_breakpoint_graph.reconstruct(args)
  File "/ngs/coral_soft/CoRAL/src/infer_breakpoint_graph.py", line 1630, in reconstruct
    b2bn.find_amplicon_intervals()
  File "/ngs/coral_soft/CoRAL/src/infer_breakpoint_graph.py", line 224, in find_amplicon_intervals
    self.find_interval_i(ai, ccid)
  File "/ngs/coral_soft/CoRAL/src/infer_breakpoint_graph.py", line 362, in find_interval_i
    if i in self.chimeric_alignments_seg[chr]:
KeyError: 'chr5'

Output from winnowmmap

m84126_240514_124343_s2/259004013/ccs   16      chr20   60748778        60     3              2149S181M1D119M1D157M1D80M1D441M1I92M1D1211M1D496M1D286M1D200M1D230M1I102M1D466M1I141M1I337M1D396M1D111M1I187M6228S    *       0       0       GCCTTGCCTGCTGGGC NM:i:20 ms:i:10346      AS:i:10346      nn:i:0  tp:A:P  cm:i:169        s1:i:2190       s2:i:0  de:f:0.0038     SA:Z:chr20,60748778,+,32229S5235M10D6151S,60,32;chr20,60748778,+,1112S5123M7D37380S,60,24;chr20,60748778,-,1034S5124M6D37457S,60,35;chr20,60748778,+,18055S4288M2I21270S,60,12;chr20,60747102,+,11464S3447M4D28704S,60,8;chr20,60747102,-,11384S3446M5D28785S,60,16;chr20,60748778,-,17974S3310M22331S,60,13;chr20,60747102,+,25639S3150M7D14826S,60,8;chr20,60747102,-,25558S3150M7D14907S,60,8;chr20,60747102,-,42501S1114M6D,60,6;chr20,60747102,+,42579S1034M3D2S,60,8;     rl:i:11102

Output pacbio
m84126_240514_124343_s2/165615159/ccs   2048    chr1    19746   60      3994S30=1X81=1X83=1X112=1X43=2D203=1X16=1X20=1X124=1X44=1X66=1X149=1X42=9542S *       0       1026    TGTTCCTTGTAAAT IDDI+D2DDI+IDDIIDII  ML:B:C,21,113,104,179,70,105,198,44,42,228,148,183,70,175,95,67,84,68,234,208,25,128,22,67,74,134,43,194,48,169,50,224,23,185,204,139,101,73,245,209,44,68,40,15,232,103,128,243,148,67,205,100,97,123,230,96,150,217,146,64,17,187,247,63,110,121,235,29,51,43,228,19,165,77,231,141,19,63,219,100,129,36,28,25,12,31,53,87,53,112,112,25,128,119,42,9,166,43,90,40,121,166,99,139,95,86,155,75,65,41,198,68,60,10,77,36,151,21,82,67,103,87,3,60,62,30,15,82,10,63,71,126,48,80,161,139,45,73,71,67,92,144,225,160,138,15,32,76,43,59,12,156,160,31,154,49,60,17,240,40,104,147,170,114,128,15,12,79,178,25,121,180,189,168,73,35,136,51,21,31,42,25,16,39,139,59,103,39,60,70,47,15,111,73,138,154,63,83,34,218,158,37,72,26,46,131,38,91,141,103,19,74,92,139,17,142,7,95,52,33,52,79,53,71,128,199,55,94,88,35   MM:Z:C+m?,14,28,33,5,7,9,9,0,16,48,7,27,5,2,27,14,0,8,10,6,2,28,18,18,15,9,15,4,15,29,5,9,27,20,28,33,5,7,9,9,0,16,10,18,15,9,16,44,13,4,3,9,27,20,28,32,5,6,9,9,0,16,48,7,27,5,2,27,14,0,8,10,6,2,29,18,18,15,9,15,4,15,29,5,9,27,20,16,15,0,8,10,6,2,27,18,18,15,9,17,44,13,4,3,9,27,20,28,33,5,7,9,9,0,16,44,0,6,25,4,32,14,0,8,10,6,2,28,18,18,15,9,15,4,15,29,5,9,27,20,27,33,5,7,9,9,0,16,10,18,15,9,16,44,13,4,3,9,26,20,28,32,5,7,9,9,0,16,44,0,6,25,4,32,14,0,8,10,6,2,28,17,18,15,9,15,4,15,29,5,9,27,20,16,15,0,8,10,6,2,28,18,18,15,9,17,44,13,4,3,9,27,20,28,32,5,7,9,9,0,16,44,0,6,25,4,32,14,0,8;    ac:B:i,4,1,0,7  bx:B:i,18,0     ec:f:5.87705    ma:i:2np:i:5  rq:f:0.998114   sn:B:f,12.4004,18.8046,4.50223,7.72926  we:i:4525584    ws:i:41102   zm:i:165615159   qs:i:18 qe:i:14578      bc:B:S,3,3      bq:i:100        cx:i:4  bl:Z:ATCACGCACACACGCGCG       bt:Z:   ql:Z:IID+7D<2<2+DDDDDDD qt:Z:   ls:B:C,137,162,98,99,164,51,45,45,51,162,98,113,100,163,99,115,97,255,164,108,101,97,100,132,164,99,53,109,99,147,5,0,0,164,113,53,109,99,147,5,2,9,162,113,108,178,73,73,68,43,55,68,60,50,60,50,43,68,68,68,68,68,68,68,162,115,113,178,65,84,67,65,67,71,67,65,67,65,67,65,67,71,67,71,67,71,166,110,101,115,116,101,100,192,164,112,53,109,99,0,163,112,115,97,0,163,115,115,97,0,165,116,114,97,105,108,192        RG:Z:a9765f35/3--3    mg:f:98.8293    NM:i:13 rm:i:1  SA:Z:chr1,190269,-,7008S112=1X83=1X112=1X43=2D479=1D157=1D11=2I22=6528S,60,3;chr1,190531,-,5188S91=2D359=1D532=8390S,60,0;chr1,190269,-,264S112=1X83=1X112=1X43=2D190=1D329=1D121=1D13=1I15=13274S,17,3;chr1,190531,+,5548S91=2D279=1D195=1I414=8032S,34,0;chr1,190269,+,622S88=1I24=1X83=1X112=1X34=1D7=1D1=1X249=1I121=1D208=1D91=12914S,7,4;chr1,190531,-,11935S59=1D31=2D62=1D12=1I814=1646S,60,0;chr1,190531,+,12297S91=2D266=1D208=1D286=1D126=1286S,60,0;chr1,190451,+,7552S14=1X112=1X43=2D40=1D627=6170S,27,2;chr1,190445,+,10919S5=1X14=1X112=1X43=2D19=1I286=1D365=2793S,25,3;chr1,190451,-,10566S14=1X112=1X43=2D106=1D93=1D470=3154S,22,2;chr1,190451,-,13938S14=1X112=1X7=1D35=2D103=1D349=,27,2;chr1,19928,+,14296S14=1X34=1D77=1X43=2D74=20S,1,2;

Why do you need to use pbmm2? Can you show the first few lines of output from pbmm2 and winnowmap?

@gringer I have updated the question and reduced the nucleotide and quality columns from pbmm2 and winnowmap output.


text:
Novogene sequencing archives
We have recently received the sequencing files from this company. For each sample, two files with extension Bam and bam.pbi are available. I would like to ask what is the meaning of the numbers behind the sequence and the quality for a specific sequence.

Please don't post images of text, and especially not gigantic images where the text is just a tiny part of it. Open the file on a computer not a phone and paste an example of its contents. I can't understand what it is from that image and you also haven't told us which of the two file types you are showing us. Please [edit] your question and clarify.

I have never heard of "PFI" files. BAM files are a format that contains the reads either with out without mapping information. What exactly you have, and how it was generated and processed you must ask the company. Only they know. Call or email them. You pay for it, so make use of that.

Most likely a pbi file? bam index?

source: Stack Exchange

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

sort: votes


text:
Tidying MEGAN Taxonomy for Python/R Analysis
I am analysing some WGS data in MEGAN and would like to do some additional analysis in Python/R
I am having trouble Tidying the Taxonomic data in a format which would be conducive to this. Originally I tried exporting Taxon Path to count and using the Split function in Pandas to break up taxonomy into new columns as this worked well for KEGG Pathways. Due to varying specificity of taxonomic assignments this didn't work.
I've had more luck Exporting Data as a BIOM file format and converting it to TSV resulting in a file with the following column containing all taxonomic info (See below)
Ideally , I'd like to generate columns for; Domain, Phylum , Class , Order, Family and Species by accessing the  features contained in Taxonomy. I'm unaware of any functions which would let me do that.
an idea of what I'm trying to achived is

Domain
Phylum
Class
Order
Family
Species

Bacteria

uncultured bacterium

Bacteria
Bacteroidetes
Bacteroidia
Bacteroidaceae
Bacteroides
Bacteroides acidifaciens

From

Taxonomy

d__Bacteria; s__uncultured bacterium

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides acidifaciens

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides caccae

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides cellulosilyticus

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides faecis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides finegoldii

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides fragilis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides intestinalis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae;

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae;

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides salyersiae

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides stercoris

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides thetaiotaomicron

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides uniformis

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides xylanisolvens

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides cellulosilyticus CAG:158

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__uncultured Bacteroides sp.

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. 2_2_4

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. 3_1_23

d__Bacteria; p__Bacteroidetes; c__Bacteroidia; o__Bacteroidales; f__Bacteroidaceae; g__Bacteroides; s__Bacteroides sp. D2

Thanks @dunc4n, pandas is a good solution. It's fairly easy. I will not have time in the next two weeks to answer this, particularly coding from an inline table. If you supplied some code that would be helpful, particularly the data you want as a Python list rather then a inline table.e

Save the following as a python script (e.g., named splitty.py). Make the script executable by doing chmod +x splitty.py. Then, if your MEGAN data are in a file called taxonfile.csv, run splitty.py by doing ./splitty.py taxonfile.csv. Get help on the script by doing splitty.py -h.
#!/usr/bin/env python3

"""
Splitty

Use this to split up a taxonomy table
"""

def argsparser():
    """
    argsparser
    """
    import argparse
    parser = argparse.ArgumentParser(prog="splitty",
                                    formatter_class=argparse.ArgumentDefaultsHelpFormatter,
                                    description="Uset splitty to get taxon table.")
    parser.add_argument("table", help="Taxon table in semicolon separated format.")
    return parser

def parseline(line_in):
    """
    Parse each line_in as return a dict
    """
    transl = {"d": "Domain",
             "p": "Phylum",
             "c": "Class",
             "o": "Order",
             "f": "Family",
             "g": "Genus",
             "s": "Species"}
    return {transl[val[0]]: val[1] for val in [i.split("__") for i in line_in]}

def main(args):
    """
    The main routine.
    """
    import pandas as pd
    df_list = []
    with open(args.table, "r") as input_handle:
        for line in input_handle.readlines():
            linelist = parseline(list(filter(None, [i.strip() for i in line.rstrip().split(";")]))[1:])
            df_list.append(pd.DataFrame(linelist, index=["Taxon"]))
    print(pd.concat(df_list).to_csv(sep="\t", index=False))

if __name__ == "__main__":
    main(argsparser().parse_args())

Thank you this got things sorted great !

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


text:
Predicting gene expression, RNA-seq, from regulatory signals, ChIP-seq
For a given cell line, I have ChIP-seq and RNA-seq data. I am trying to construct a model for predicting gene expression from ChIP-seq. Basic  models  I have seen [a very simple model is described in  https://www.pnas.org/doi/10.1073/pnas.0904863106 ] generally rely on chip-peak-to-gene-distance:
Gene_Expression(gene)=F_gene(ChIP-Type,tag_count,distance_to_gene)
This peak to gene distance dependence seems over-simplified, so I am looking for alternative, possibly more nuanced, quantitative models to describe ChiP (TF) to gene expression dependence.

I predict that this question will receive many "needs more focus" votes. The entire problem of peak to gene association is very unsolved. There is reviews that summarize the different approaches that have been developed over the last two decades. Please read them first, and then ask a precise technical question.

Thank you @ATpoint, I recognize  the question is posed generally - and that is somewhat the point here - that you've touched on with your "very unsolved" characterization. What I have seen, in the few articles I have read, are arbitrary specialization of  models without satisfactory justification - getting the sense of a struggle to formalize the variations of functional connection types between Gene expression and regulatory elements/signals. For this reason I am looking for broader/general input.


text:
How to identify genomic regions / peaks associated with enhancers (TF binding sites)? Is there a tool or a formal recipe?
I am attempting to identify DNA sequences/regions associated with enhancers in my NGS sequencing data. How can I identify (from peak files/counts for chip-seq, atac etc.) genomics regions/peaks associated with enhancers? I know proximity to genes and possibly histone (i.e. h3k27ac) activity near known genes is possibly an indicator, but is this the only way and how is the computational process formalized/captured? What are the specific set of tools involved? i.e. Peak annotations etc.? Thank you.
As I understand it in my limited bio knowledge, enhancers are at varying distances from genes and sometimes they are multiples (super-enhancers). Are all enhancers already "known/described" (akin to genes)? What is/was the process of finding them (human manual procedure)? It's that more manual/heuristic process (without ML) I am looking for so we can consider formalizing and automating it.
Based on a very insightful answer: I think some input/clarity on the following bit would make the answer(s) even more valuable:
Possible list of steps towards this include following tools/functions:

Bowtie (to align fastq to a genome, mouse, human etc.)

Homer/makeTagDirectory (Sep files, associate chr/coords with tags/sequence reads)

Homer/findPeaks (identify and list genomic regions/sequences with count values above a threshold)

Homer/pos2bed (convert peaks to bed file, genome browser compatible form)

Homer/annotatePeaks.pl (Runs on *bed file - Don't know much about this, but intuitively associate a label with peaks found in 4?)

Question: Is there a step 6? Where do we go after annotatedPeaks for example to get to enhancers?

If you have reads from chip-seq or atac you can use tools like homer or macs to identify peaks in your data. Once you have the peaks, these programs allow you to do motif discovery to annotate your peaks with e.g. what transcription factor binding sites they contain.

Thank you @Niek. Very insightful contribution, that has allowed me to add specificity to the question which I have updated with your input.


text:
How to chop fasta / bed /peak file genomic segments into smaller fixed or custom genomic intervals
I am trying to split my .bed/.fa file genomic segments into arbitrary smaller overlapping intervals:
Consider a typical line/row in my *.bed file as follows:
chrom. $\ \ \ \ \ $  start $\ \ \ \ \ \ \ \ \ \ \ \ $ end
chr7 $\ \ \ $   127471196 $\ \ \ $    127472363$\ \ \ $      ...
The start and end genomic positions include 1167 base pairs for this line. I need to extend this line, for example, into multiple lines where each line covers exactly 500bp like below:
chr7 $\ \ \ $    127471196 $\ \ \ $    127471696 $\ \ \ $       ...
chr7 $\ \ \ $    127471863  $\ \ \ $   127472363 $\ \ \ $       ...

While each line covers 500pb, between the two new lines they both include the same "range" of genomic positions as the original line including some overlap which is fine (and desired)
Is there an existing tool (perhaps part of bedtools?) that can do something like this?
I know there is the bedops tool that does ALMOST this using the --chop option, where it takes the original segment x1:x2 and chops the original segment into specified fixed bp length intervals, however, if there is any left over towards the end is creates a new non-standard length segment (i.e. x1:x2=10bp and chop length=3, it will give you 4 segments that are length 3, 3, 3 and 1 - what I am looking for is something that adjusts for the last remainder and outputs segments of length 3,3,3 and 3 where the last length 3 segment is right adjusted, while the first three segments are left adjusted or better yet it gives the user the choice of precise redefinition of segment coordinates y1,y2,y3 and y4).  Thank you.

You can run --chop with -x to excise the last piece, adding the --stagger to size the chopped pieces as needed, e.g.:
$ echo -e 'chr1\t0\t10' | bedops --chop 3 --stagger 2 -x -
chr1    0       3
chr1    2       5
chr1    4       7
chr1    6       9

Adjust parameters to taste.
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#chop-w-chop
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#stagger-stagger
https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#exclude-x

Thank you @alex - the specific issue with bedops --stagger option is that I have not been able to feed it a negative number, that's what would allow for "overlaps" which is critical!

As shown in the example above, --stagger can give you overlapping elements. However, you also can pipe results to bedops --everything --range N to extend edges by N bases, e.g. ... | bedops --chop X -x - | bedops --everything --range Y - > answer.bed. You'd probably need to use coordinates from fetchChromSizes or similar to trim edge-of-chromosome windows, but that is easy to do via piping to awk or the like.


text:
STAR vs Bowtie2
What is the fundamental difference between STAR and Bowtie(2). Specifically, what is the difference in their final output (regardless of run-time differences, speed, memory usage etc.). Both seem to work on the same kind of input *.fastq files without discriminating what they are associated with and generate an output version of the *.fastq files. Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used, and vice versa? Thank you.

What is the fundamental difference between STAR and Bowtie(2).

STAR for mapping spliced (i.e. with introns) short RNA-seq reads against a genome. Bowtie2 for mapping short reads without splicing.

Would my downstream analysis be incorrect if I use Bowtie2 where STAR was supposed to be used (?)

Yes. If you use bowtie2 to map spliced short reads against a genome, STAR-based downstream analysis is likely to be wrong.

... and vice versa?

If you map genomic reads with STAR, the result will be suboptimal. Use the right tool for the right task.

The seeds are clustered together by proximity to a selected set of ‘anchor’ seeds in the 2nd phase of the STAR algorithm,  but this •clustering/stitching/scoring process• is not present in Bowtie2.

mapping by performing a two-step process


text:
PWM (4xM matrices) to Sequence Logo Visualizations in MATLAB or Python
Is there a MATLAB, or else Python (not R) tool for visualizing sequence Logos from probability weight matrices (PWMs)?
This seems like a basic and simple tool to implement, yet just above the threshold in the space of tedious to write my own from scratch - I was surprised I couldn't find a basic function like a
visualizePWM(pwm_matrix)

that would take a pwm_matrix that is a (4 by m) matrix and spit out the visualization of an m-length sequence motif/logo like:

I see all kinds of over-complicated bloated tools that I thought I ask in case I am missing something.

answer from @M__, converted from comments:
This has been answered on stackoverflow using matplotlib for seq logo plot:
https://stackoverflow.com/questions/42615527/sequence-logos-in-matplotlib-aligning-xticks
IMO its only useful for presentation, pretty pictures, thus there's no need for someone to write a library (which is what you're asking). If I'm using a mutation matrix (which I do), I want the raw numbers and significance values I don't want pretty graphs getting in the way.
Conceptually you would move from 2D to 1D, without subplots (which is doable), i.e. its a matrix. My view on code availability is not universal. I've had requests for separate CLIs (useful for someone clearly) ... which just never happened.

The solution proposed in the link is focused primarily on a tool that is captured as the pyseqlogo packet at the end. I have spent a good part of the day trying install the packet and get it to work, but to no avail it's failed with too many errors and complications to list :(.

After much searching and trial and error, located possibly the optimal tool for my question:
seqlogo
https://github.com/betteridiot/seqlogo
Surprisingly not only seqlogo is an elegant and minimal solution in Python (and R) I have discovered it has a relative implementation in MATLAB! It installed and worked the first try - beautiful output with all the output options I needed - rock solid.
Python one liner (not counting the example pwm input generation)
random_ppm = np.random.dirichlet(np.ones(4), size=16)
ppm = seqlogo.Ppm(random_ppm)
seqlogo.seqlogo(ppm  ,format='png',filename='out.png')
Produces the following ('out.png')

And here is MATLAB Port:
https://www.mathworks.com/help/bioinfo/ref/seqlogo.html


text:
Common tools for generating results and visualizations for bioinformatics and genomics publications
I am a computational scientist with background in computer science and engineering however very new to the field of bioinformatics and computational genomics. Looking at papers in nature, PLOS, Cell, PNAS etc. within the topic of computational genomics, I see a very similar set of very sophisticated graphical results, ranging from clustered/sorted heat map representation of gene expressions to log2 fold changes, volcano plots, circular graph network visualizations etc. Coming from an engineering background I would often generate customized tools for visualizing results myself as part of each paper, albeit not at the sophistication level of results I see in bioinformatics. But here, it seems there is a set of tools standardized where many are using them to display their results. Is there a summary of these result generation/visualization tools along with the references to their implementation/usage (opensource I assume)? Something like that could be very helpful to newcomers like myself. Thanks in advance.
Below is an example /elaboration based to guide the question further into specifics:
Consider the example paper:
https://doi.org/10.1016/j.cell.2018.04.018
This paper contains a principle set of results/visualizations I see almost universally in many nature/cell computational genomics papers. I am wondering if there exists an integrated library either in python, R (or Matlab for example when they create toolboxes and integrate a commonly used subset of tools in an area) that everyone taps into for generating these type of results. The "principal set" that I am referring to include:

Genome browser track snippets: Figure 1A
log2Fold x-y scatter plots : Figure 1B
Hierarchically clustered/sorted gene exp heatmaps: Figure 1C,1F
Venn diagrams : Figure 4F
Hi-C type square ligation heat-maps: Figure 6A

Figure 1A: These genome browser snippets are either plain screenshots from applications like IGV or created from within R with software like gviz.
Figure 1B: This looks like base R plotting routine. Like plot(x,y) and then some coloring and legending. No particular "tool" here other than really just base R or ggplot2.
Figure 1C/F: Heatmap packages include ComplexHeatmap and pheatmap. I personally use the former as it is versatile with extensive documentation.
Figure 4F: Venn diagrams, just google "Venn" and "R" for several packages that do it in R. Same for other languages. There are plenty of packages. No preference.
Figure 6A: These plots are often returned by specialized HiC software such as HiCexplorer which is powerful with lots of documentation.


text:
How much of the space of motifs have been discovered?
I notice when using tools like HOMER there are two different types of outputs: de novo motifs vs known motifs and  it seems no matter what the data is (type of ChIP etc.) there is always a "de novo" motif set in the results.
I was hoping to reach out to experts on insight about how much of the space of motif based structure, within a genome, has been identified. Is there any consensus on what portion of the motif like patterns, for example in human or mouse genome, have been identified? When I find de novo motifs using HOMER in my ChIP-seq or ATAC peaks, what are the chances that in my common-place ChIP experiment (PU.1 for example) I find de novo motifs, not recorded in the databases (i.e. JASPER).

Thanks. I previously closed this question mostly due to the "why" in the title. I've reopened the question because it looks more likely to encourage specific answers.

De novo means that the software infers presence of motifs based on the provided sets of DNA sequences, most commonly peaks from ChIP-seq and similar experiments.
So essentially it stacks the sequences on top of each other and checks for presence of conserved nucleotides that form sequences of certain lengths, for example < 10bp.
Known means that each peak is checked whether it contains a known motif (so is similar to a given known motif) and then in the end a per-motif statistic is calculated to decide whether enriched or not.
If the investigated protein is well-studied then de novo and known is likely similar or identical. De novo does not mean that found motif is novel in terms or unknown or undescribed. It just means that it is derived from the data directly rather than matching known motifs onto the peaks. It's simply two different ways to approach the question. Hope that makes sense.

Thanks @Atpoint, greatly useful and helpful answer.  Part of my question remains, albeit in a HOMER independent context, as to how much of the motif space remains undiscovered, if any?

Probably quite some. Imagine just for how many transcription factors there are no antibodies available.


text:
How do you convert Raw Alignment Score to Bit Score?
I'm coding a pipeline where I make a lot of pairwise alignments, and I end up with raw alignment scores. But, I really need to look at my results in terms of bit scores.
I know that the formula is:
$′ = (ƛ* − ()) / (2)$
But, I don't know what ƛ or K would be. I'm currently using BLOSUM62.
I've looked around online for how to look up (and even empirically deduce) lambda and k, and I all I've come across were convoluted explanations.
I also know that this question from almost 3 years ago attempts to address the same topic, but the best response only leads to an NCBI link that requires login credentials.
Any help would be much appreciated.

Answer from @terdon, converted from comment:
The link just points to the NCBI ftp site where you can download the source package and look at the code. There is no need for login credentials, you can log in as anonymous. It's just that the version has changed so that answer's URL no longer works. Try this instead:
wget https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST/ncbi-blast-2.11.0+-src.tar.gz

Then, extract the files and look at the file c++/src/algo/blast/core/blast_stat.c which, again according to the answer you linked to, should have the lamda and k tables at line 589.


text:
Simultaneously get data from multiple applied gates in flowCore
Using the Bioconductor flowCore package, I'm applying two parallel and non-overlapping gates to a gatingSet directly under "root":
library(flowCore)

# Import data
file.name <-system.file("extdata","0877408774.B08",package="flowCore")
x <-read.FCS(file.name, transformation=FALSE)

# Create flowSet
fs <-as(list(x,x), "flowSet")

# Gate
library(flowWorkspace)
gs <- GatingSet(fs)

gate1 <- rectangleGate(filterId = "Gate1", list("FSC-H" = c(300,600), "SSC-H" = c(50,300)))
gs_pop_add(gs, gate1, parent = "root")

gate2 <- rectangleGate(filterId = "Gate2", list("FSC-H" = c(700,800), "SSC-H" = c(500,600)))
gs_pop_add(gs, gate2, parent = "root")

recompute(gs)

gs_get_pop_paths(gs)
# [1] "root" "/Gate1" "/Gate2"

# Plot
library(ggcyto)
autoplot(gs, gs_get_leaf_nodes(gs))

# Extract gated data
gs_pop_get_data(gs, "Gate1") # Works
gs_pop_get_data(gs, "Gate2") # Works
gs_pop_get_data(gs, c("Gate1", "Gate2")) # Does not work
gs_pop_get_data(gs, "Gate1", "Gate2") # Does not work

I want to get the data from the union of gates 1 and 2 but I can't figure out how. The flowCore documentation touches on this topic, but the examples shown are for filters, and I haven't figured out how to get it to work with gates.
If there is a more straightforward way of doing this using another package that would also be of interest, as long as it is compatible with flowSet/flowCore objects.


text:
Calculating the segregating sites via tskit versus counting
Can you help me understand why I am getting different values when I calculate the segregating sites from tskit api versus manually counting? Thanks!
Example data generation:
import sys, tskit, msprime

ts = msprime.sim_ancestry(samples=3, ploidy=2, sequence_length=10, random_seed=2)
ts = msprime.sim_mutations(ts, rate=0.1, random_seed=2)

Example results: I get 5 from API and 6 from manual counting.
>>> ts.segregating_sites(span_normalise=False)
5.0
>>> n=0
>>> for tree in ts.trees():
...     for site in tree.sites():
...         for mutation in site.mutations:
...             node = ts.node(mutation.node)
...             pop_id = node.population
...             n+=1
...
>>> n
6


text:
Nextflow: How to initialize channel from CSV file but only after CSV file created?
Process A runs a python script which generates a CSV file in some absolute file path say "/mnt/x/blah.csv".
Process B needs to consume the CSV file line by line, after process A has finished.
The trouble is that the channel to input to B is initialized at the beginning before "/mnt/x/blah.csv" is created by A so it fails. This occurs even if I collect the result of process A as input to process B (because channel is initialized regardless)
nextflow.enable.dsl=2

def get_the_csv_records() {
    Channel
        .fromPath("/mnt/x/blah.csv")
        .splitCsv(header:true, quote: '\"')
        .map{ it.TheFieldIWant }
}

process process_a() {

    output:
    val("I'm done!")

    '''
    python.exe ./make_the_csv_file.py --path=/mnt/x/blah.csv
    '''
}

process process_b() {

    input:
    val(previous_step_is_done_dummy_variable)
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a()
    process_b(process_a.out.collect(), get_the_csv_records())
}

This fails immediately with "No such file: /mnt/x/blah.csv" because the channel is initialized before process_a has even started.
Same behavior occurs when using concat:
    process_b(process_a.out.collect().concat(get_the_csv_records()))

Question:

How can I make process_b() only read the input channel from the CSV file after it is created?

Alternatively, can I make process_b() trigger when process_a() is finished some other way?

Alternatively, can I add an intermediate step which reads process_a.out.collect(), then reads /mnt/x/blah.csv and outputs it into some existing channel which process_b() takes as input?

Another approach - what is the way to create a channel which outputs the results of splitCsv, when the file is not known at compile time?

Assuming you are using latest version of Nextflow which support the DSL 2 language, and assuming your python process output the file blah.csv, then the easiest way to move forward will be:
process process_a() {
    output:
    path("blah.csv")

    '''
    python.exe ./make_the_csv_file.py --path=blah.csv
    '''
}

process process_b() {

    input:
    val(the_field)

    '''
    echo ${the_field}
    '''

}

workflow do_a_then_b{
    process_a() \
        | splitCsv \
        | map{ it -> it.TheFieldIWant } \
        | process_b
}

You can also refer to the nextflow google group.


text:
seqff script got trouble after samtools treatment
I'm trying to study fetal fraction by using a seqff script from here. The instruction says that I need a headless sam file to work so after alignment with bwa, I process my sam file by using samtools:
samtools view -S -q 20 -F 4 target.sam > treated.sam

This process also remove the @SG header since there is no reference file. I also use grep to remove header of the target.sam got after alignment for comparison. However, when I compare, I noticed that ENet calculation on the treated.sam returned with NA results, so half of the method to calculate fetal fraction of this script was gone. I have tested it in multiple cases and it seems that without using samtools, the script worked well but with samtools Enet calculaton not worked. Can someone tell me what is wrong with the script? I'm using R 4.0 and all modification I made to the script is on the args since it cannot recognize args I input with anymore.
p/s: this is few lines of headless target.sam made by using grep
W2003120038L1C001R0050000009    0   chr10   6490952 0   50M *   0   0   GAGCAGAAAAACTGGAAACTCTAAAAAGCAGAGTGCCTCTCCTCCTCCAA  IBGIBAB<ABBCA9IB@BHDIA=ABBB8IB9A;DGII@I8IIAII<IDCA  XT:A:R  NM:i:0  X0:i:151    XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26

And this is from treated.sam, the sam file be processed by samtools
W2003120038L1C001R0050000037    0   chr5    172722817   37  50M *   0   0   GTCCTAAGACAAAGACCTCAGTGCTTATAATGTAAGTGGTCACTGGGCTG  ICIICBB:BIBBBIAICDIAIAIICDBDAABIDBBICEIAIAFDH?IIAI  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7T42
W2003120038L1C001R0050000062    0   chr5    84724587    37  50M *   0   0   TAAGAATTGTGCATAACCAGTAGCCCATGTTATCTTGCAGAAGCATTTGT  ABBI@BDDIDIIADABHEAIC@IIIBBAID@ADIDCA:BC=BIIAAAAIC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000067    16  chr5    45648067    37  50M *   0   0   TTTATTTAAATACAAGCCAATTAAATGATTTTATTTTACATTTTCTAACC  AABAAABCAAABH@CFHIDDB@BCDAID@BBACBB@BBIDB@2BIBC@II  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000184    16  chr4    103629659   37  50M *   0   0   CAGTGCTAAGAGGGAAGTTTATAGCAACAGATACTTACATCCAAAAAGTA  FDFBIIADAFDICICCIAABDBAFIADIDFDABIBA@I@CI3?DCCDIBC  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000211    0   chr3    157734199   37  50M *   0   0   TCTGAAATTGATCTTTATCCCACATCCACGTATCCTAAAATAATAAAACA  AFDBBBB*CIBDFACABDIIBB<1CI=BFIDBAFEABBAADBBABAABIB  XT:A:U  NM:i:2  X0:i:1  X1:i:0  XM:i:2  XO:i:0  XG:i:0  MD:Z:7G15T26
W2003120038L1C001R0050000212    0   chr4    157022243   37  50M *   0   0   AAAAATAAAAAAGCATTCTGCCCAGATCTCTCAACATTATGTTATTTCTT  4<AA?D2,AAA@IF@DDIDIIDIA2AD=DBAF@BIADA<DICD@DADFD@  XT:A:U  NM:i:0  X0:i:1  X1:i:0  XM:i:0  XO:i:0  XG:i:0  MD:Z:50
W2003120038L1C001R0050000258    0   chr13   106093727   37  50M *   0   0   AATGGTTACTTTATGTTTATATTTCAACAGGGAATATTTGAAATATGGCA  97?IIDA*FDAD<@H@AD<@<?AAG8??,I;I@:DADDAG7@@7(;IH6?  XT:A:U  NM:i:1  X0:i:1  X1:i:0  XM:i:1  XO:i:0  XG:i:0  MD:Z:7C42

@gringer I have added few lines of both file on my post. Also I have to perform quality filter to filter out bad alignments results

If the grep one works, and the non-grep one doesn't, my first step would be to exclude any confounding factors (like quality filtering) to try to get them both to look the same (which they don't, at the moment). Is there any improvement if you don't include the quality filter (either in the files looking similar, or the script working better)?

If you think there is a formatting problem, you might need to put small subsets of the samtools file through until you can identify the problematic line.

@swbarnes2 thanks, do you know any command or program that can check for the different? Because the each sam file is more than 4GB, finding what is different by compare with text editor is impossible to say at least

One you figure out the problem in one file, you should be able to fix it for the others.


text:
How to calculate RMSD and binding energy from the already docked ligand-receptor
I have the following PDB file that looks like this.
It is downloadable here. This PDB already contains ligand and receptor together.
The red molecule is the ligand and the yellow is the receptor.
They're already docked with ClusPro.
My question is how to calculate the RMSD and binding energy from it?
Any existing tools (with or without command line) that can help?
I don't have any data to try it out at hand but have a look at the documentation I linked. There you can find ways of constructing an Ensemble which can then be superposed and RMSDs calculated, even pairwise.
Reference is the red one (receptor). I have hundreds of other complexed ligands. I want to calculate how good are those complexes. What's the reasonable measure to do that?
Example code using Prody? please :-)

Could you figure it out how to calculate the binding energy? which software did you use? Thanks,

NB. This is a long comment as opposed to an answer which requires a colossal working chunk of code.

Control
Science is all about controls failing. A great control for docking is the barnase-barnstar complex (PDB:1BRS): both small, but form a strong interactions. If this is done, the RMSD against the original structure is very meaningful and useful.
RMSF
A RMSD is a metric of how much two structures differ. In an ensemble, the metric of how much a structure and the ensemble average is called a RMSF (root-mean square fluctuation). Here the RMSF is sought as the RMSD from first/top pose or worse the manually thrown together pre-docking pose does not mean much.
Clusters
The decoys (PBDs trails) need to be clustered into different ensembles first, before a RMSF can be calculated.
∆∆G and water
For the calculation of difference in Gibbs free energy between the docked minus the docked poses, the structures need to be energy minimised with the forcefield scorefunction that is to be used for scoring.
Water is very important for protein-protein interactions, but adds a lot of complication. Say Gromacs is used, gmx energy will return the single state (conformer) Gibbs free energy, but it requires a water box. For the difference, the water box of undocked needs to be the same size and preferably with the same amount of molecules. A much simpler solution is using an implicit solvent model, e.g. Rosetta. As a bonus, the energy minimisation (relax) is relative fast. And specific explicit waters can be added. However, calculating ∆∆G for every decoy will take hundreds of CPU-hours and is a bit nonsensical.
A better approach is to look at residue conservation (consurf), cross-linking mass spec data if available, and if human post-translational modifications (phosphosite) and healthy variants (gnomAD). But ultimately wet lab validation by ITC of different mutants/labelling based upon the docking results etc. is required.


text:
How to make a video of a gromacs gmx dump standard output with VMD
I did gmx trjconv on GROMACS trr file to generate a .xtc file, and then used gmx dump and saved its standard output in a file. Then I tried to make a movie with VMD using the file. As a result, it does not work properly. The error in the terminal is as follows.
gromacsplugin) Cannot read header fromm '/Users/kondouakira/mount/Lab/1ns_sim/md_0_1_noPBC3.trj', no error

I use following command.
gmx trjconv -s md.tpr -f md.trr -o md_0.trr -pbc mol -ur compact -b 1 -force
gmx dump -f md_0.trr > md_0.trj

No need to use gmx dump, you can simply load a xtc or trr trajectory in VMD using:
vmd path/to/gro_file.gro other/path/to/trajectory_file.xtc

As always you can also use the GUI to first open your gro file, then right click on the molecule and select Load data into molecule before selecting your xtc or trr file.
After putting the representations as you see fit, if you want to export this as a movie (GIF, MP4, others...) you can go in Extensions -> Visualization -> Movie Maker. See Trajectories and Movie Making VMD official tutorial for more details.

Thank you for your answering! But I want to draw a video using a trj file.There is also a trj file option in the VMD file format. Therefore, I would like to know how to perform this to create a video.

I don't think that the output of gmx dump can be read by VMD: the goal of this command is to provide a human-readable output, not to convert one trajectory format into another.
trj files are NOT the output of gmx dump. trj files were used in "old" versions of GROMACS to store trajectories (see trj file format documentation of GROMACS 5), and this is what VMD is awaiting, not the output of gmx dump. This is also why it tells you the header is wrong.
In recent versions of GROMACS trj files are no longer supported, probably as they were not in a binary format, hence took too much disk space (see Trajectory files formats in GROMACS 2021 documentation). Converting a .trr in a .trj was done by gmx trjconv, but I don't think this is possible anymore (see gmx trjconv supported trajectory formats in GROMACS 5 vs GROMACS 2021).
Why would you absolutely need gmx dump output to be read by VMD? The only reason I can see is a handmade script doing modifications to your trajectory. In this case I would recommend to instead look into some molecular dynamics python libraries, like MDAnalysis, PyTraj or MDTraj, which will also allow you to make minute modifications to the trajectory file while still reading it in its original format. You can also look harder into GROMACS tools, as they are quite complete.


text:
How to show that one docking is better than another with Pymol visually
I have the following docking results:
Dock_1 has better energy score than Dock_2.
What I want to do is to show with Pymol that Dock_1 is better visually.
Here I use hydrogen bond, highlighted with Yellow (Dock_1) and Gray (Dock_2).
But as you can see, with that I can't show that Dock_1 is clearly better.
That's why I'm looking for better alternative than hydrogen bond.

In the last paragraph of this post are some suggested things to check that may be useful: https://bioinformatics.stackexchange.com/questions/16192/how-to-calculate-rmsd-and-binding-energy-from-the-already-docked-ligand-receptor

maybe you can use the "per-residue interaction scores" (wich is available via Schrödinger).
On the other hand,to show if a Docking Pose is really "better" you could show specific interactions rather than the complete big molecule.. e.g. an important salt bridge is missing. (hydrogen bonds are not that important, due to the fact that they can be bridged with water molecules)
Sorry, this is not really a solution, but I am not allowed to comment jet -.-

Totally. To further this, a salt bridge contributes something like -2.5 kcal/mol, while a h-bond -1 kcal/mol. pi-pi bond or S-pi about -1 kcal/mol, but that does not account for the improved solvatation term, which plays a role.  Coordinated waters play a massive part in many protein, e.g. G-protein alpha–beta interaction, ribosome, barnase-barnstar etc., and a water can contribute some -3 kcal/mol.

@MatteoFerla my key concern is to show that visually.  Is there  away I can show salt-bridge, pi-pi bond, S-pi bond with Pymol to show the difference of two docking?

@Lucas thanks. Is there a way to show "per-residue interaction scores" visually in Pymol? Is there away we can systematically identify specific interactions that matters? Any reference on this? I appreciate your reply.

@MatteoFerla Also, if I want to show salt-bridge with Pymol. What's the reasonable residue to choose from each ligand and receptor?

Your E looks like a ∆∆G in kJ/mol, which means that they were docked with a forcefield method. You may have among the files a per-residue split of the forces.
For PyMOL, select salty, ((byres (resn ASP+GLU around 3)) and resn LYS+ARG) or ((byres (resn LYS+ARG around 3)) and resn ASP+GLU) should work for salt bridges btween canonical amino acids.


text:
How to Read SCF file in Python?
Is there any way that I can read SCF file in python like in R using sangerseqR, I have tried with Biopython, it seems it does not support this format.

You could try converting it to ABI format as asked in this question and then import that

Alternatively, bioconvert supports conversion of SCF to FASTA and FASTQ formats.

Conversion to SCF to FASTA and FASTQ formats will not help me since I need the chromatography data not just the seq

It was done using sangerseqR with rpy2 in python
    import rpy2.robjects as r
    from rpy2.robjects.packages import importr
    utils = importr('utils')
    utils.install_packages('sangerseqR', repos="https://git.bioconductor.org/packages/sangerseqR")
    utils.chooseBioCmirror(ind=1) # select the first mirror in the list
    utils.install_packages('sangerseqR')

    sangerseqR = importr('sangerseqR')
    Trace = sangerseqR.readsangerseq( file="1I1_F_P1815443_047.scf") #scf sequence file is now loaded

sangerseqR Documentation(pdf)


text:
Recode bgen files against reference genome
I have a set of .bgen files split by chromosomes, including both SNPs and INDELS (named with both rsids and chr_pos_a1_a2). I know a large portion of them are encoded on the wrong strand, and I want to recode them against a reference genome (in my case, 1KGP).
How can I recode the variants in my starting files against the reference genome?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.


text:
Precisely clipping bam file to bed coordinates
I have a bed file of coordinates and would like to clip anything outside these regions in my aligned bam files. However, trying programs such as bedtools intersect and samtools ampliconclip (on complementary bed file), they always let some reads slip through outside the ranges. I tried using hard-clip and both-ends options for ampliconclip aswell.
Example (top is clipped bam, middle is unclipped bam and bottom is bed file with regions I would like to keep):

Although most reads outside the designed amplicons are removed, there are still reads outside the regions in the Designed_pool1 file.
The problem is that the bed file contains amplicons and outside the amplicons are primers where I do not want to call variants. I cannot get the exact primer sequences because I don't know their length and the manufacturer does not state them, so I cannot use programs that require primer sequences.
Does anyone have some experience with this and know why the programs behave this way?
EDIT:
The exact commands I have used are:
bedtools intersect -a $BamFile -b $BedFile > $BamOut

and
samtools ampliconclip --hard-clip --both-ends -b $ComplBedFile $BamFile > $BamOut

as well as without --hard-clip and --both-ends options. All methods let some sequences stay out of the regions in bed files.
I am not completely familiar with the bedpe format, but I am not using paired end data. I cannot get primer start and stop positions, only amplicon positions.

Maybe you need to sort the bam file before trying to clip.  Did you try that already?

Yes, the bam files are sorted

For single-end reads, generically:
$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --fraction-ref 1 reads.bed roi.bed > answer.bed

The BED file answer.bed contains filtered reads which entirely overlap intervals in roi.bed. It contains fields required for backconversion to SAM/BAM.

I am curious whether your pileup visualization tool is respecting the clipping. It is possible that there is a subtlety of samtools ampliconclip that your visualization is not understanding, with regard to how it treats e.g. supplementary alignments or something.
Perhaps more importantly, I notice that the reads that are not being properly clipped in your visualization tend to have one end anchored in each of your BED regions.
Looking at the samtools ampliconclip docs, it seems to be assuming that reads will only be clipped from the ends. They will not be clipped if the region outside the BED is in the middle, if I read it correctly. I think that you may need to include an additional step in which you either trim reads (so that they are not long enough to bridge your regions) or you somehow postprocess your alignments to explicitly clip in the middle.
If true, this seems like a significant enough oversight of samtools that it may be worth logging an issue with the samtools maintenance team.
Also, not clear if it's relevant but reading the samtools ampliconclip docs:

adjustments to the left most mapping position (POS) will mean that coordinate sorted files will need resorting. In such cases the sorting order in the header is set to unknown. Clipping of reads results in template length (TLEN) being incorrect. This can be corrected by samtools fixmates. Any MD and NM aux tags will also be incorrect*, which can be fixed by samtools calmd.

You can convert a SAM stream (mandatorily with a header) into FASTA (multiFASTA file padded with dashes) using gofasta:
samtools view -h some.bam | gofasta sam toMultiAlign --start 1171 --end 1177 > region.aln
Beware gofasta sam toMultiAlign discards INSertions from the sample relative to the reference, so the FASTA sequences may not be exactly sub-strings of the original data. You can detect DELetions in the sample reads when - are present relative to reference.
That is not an issue with gofasta sam toPairAlign but is terribly slow and it dumps zillions of single-entry FASTA files into your current working directory. Probably would be worth to add some flags to samtools to omit unmapped reads from the input to feed gofasta with less input.
samtools view -h some.bam | gofasta sam toPairAlign --reference someref.fna --start 1171 --end 1177 > region.aln

I wrote https://jvarkit.readthedocs.io/en/latest/PcrClipReads/
java -jar dist/jvarkit.jar pcrclipreads -B test.bed  samtools-0.1.19/examples/ex1.bam


text:
Getting exact mismatch data from multiple sequence alignment
TL;DR Can anyone point me in the direction of a tool that can be used to input different viral substrain reference sequences (fasta) as well as sequence data (fastq/bam/vcf) and find the most likely substrain of the sequenced data?
I am doing multiple sequence alignment of different viral subtype references and would like to get the exact positions and base each time there is a mismatch. I have tried using different tools. For example Clustal Omega can almost give me what I seek, but I have trouble getting the exact position and nucleotide for the mismatches (the full thing is ~8000bp long):
ACTACAATAATCCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGTGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAATTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG
ACTACAATAATTCATGTATAAAACTAAGGGCGTAACCGAAATCGGTTGAACCGAAACCGG

*********** *********** ****** *****************************

E.g I would like something like:
Pos 12:
C
T
T
T
T
T
T
T
T
T
T

For each position with a mismatch.
A problem is that the sequences are slightly different lengths so something like a variant call format on the "main" strain might work, but I'm not sure how I can do this on multiple files of this length.
The end goal is actually to be able to call which substrain is present from NGS data. I imagine using variant calling on the "Main" strain, but getting from there to calling which substrain is in the data is a bit of a long way for me. If someone have any suggestions on tools for this, please enlighten me :)
Thank you very much!

I found a solution.
For variant calling on multple sequence alignment I can recommend snp-sites.
For anyone facing a similar problem of wanting to find substrains, here is an overview of useful tools:
"Computational Methods for Strain-Level Microbial Detection in Colony and Metagenome Sequencing Data"
(https://www.frontiersin.org/articles/10.3389/fmicb.2020.01925/full)

I have partly answered this for know! For anyone interested i use kalign to align all 10 strains to a fasta format, then I use snp-sites to output a vcf file that can show each mismatch in accordance to the first strain in the fasta. New problem is that the alignment "stretches" the length of this first strain, because some of the other strains are a bit longer. This means that the position of variants are with reference to the longest strain of the 10 and not the strain I want. It can work for now, but any solution to this, would be helpful.


text:
Is it possible to correct bond order and conectivity problems in an sdf file according to mol2 files?
I was doing a Virtual Screening Search using Pharmer software. The .sdf output file contains a column in mol format and another column with Root Mean Square Deviations. I can see some bond order and connectivity problems in the mol column. I do have a mol2 file that contains the original molecules and I was wondering if I can use these mol2 files, to fix this. Can rdkit work in this case?

RDKit has a one line command to fix the bond order (docs)
 mol = Chem.MolFromMolFile('wonky.mol')
 template = Chem.MolFromMol2File('template.mol2')
 newMol = AllChem.AssignBondOrdersFromTemplate(template, mol)

The problem is that the two need to be identical. I see a carboxylate in your screenshots, so you most likely did the correct thing and docked with the appropriate protonation for pH 7, i.e. as the conjugate base with a negative formal charge. If one does not then RDKit will fail.
Also, RDKit is "not great" with mol2 files. If you have the original SMILES that is better. If the problem is atom names, using a PDB is a better call. If you want to assign atom names, this blog post of mine may help.


text:
Sequence alignment with MUSCLE
For protein multiple sequence alignment (MSA), I am using MUSCLE (v5.1). The default setting without any specifications gives a 'bad' output even though my proteins are very similar (UniProt sequence and the corresponding PDB SEQRES sequences). 'Bad' is defined as unnecessary gaps are introduced into the sequences.
However when the super5 algorithm is used the MSA output is correct, which is described here https://drive5.com/muscle5/manual/cmd_super5.html
In addition, when I increase the gap penalty this issue can be solved but as far as I know, super4 v5.1 doesn't give you the opportunity to adjust gap penalty.
Documentation describes super5 quite similarly than the default algorithm, it "just" reduces runtime for long sequences.
Can someone please give me a hint what am I missing and/or why super5 outperforms default in terms of accuracy in such "easy" cases?

Hi there, would you mind making your example reproducible? Or at least show us the two alignments? I wonder how it is possible the gap penalty is not adjustable. Is there a reason why you want to use MUSCLE in particular?

In general I would agree with your comment @KamilSJaron . This particular question is an exception  because the underlying algorithm in context is very new and is precisely focused on the observation of the OP. I would not have been aware of the release of Muscle 5 without the OPs question.

Hi!
Thanks for both of you!
Now I get it. To be honest, there's no particular reason for using MUSCLE. I just found it easier to implement it in Python with many cases and the output also seemed easier to parse than lets say BLAST.

HMM Super5 is a hidden Markov model (HMM) where  gap penalties and the substitution matrix are the parameters at the centre of the model. That is the answer, thats why you can't change it - it does it itself and thats its purpose. Thus one algorithm (super5) is HMM for the gap-penalty and the other one I suspect isn't.

Explanation To explain that further HMM uses observed outcomes to estimate/optimise the hidden parameters, the gap penalty being one of the two key hidden parameters. Thus the algorithm is exploring a parameter landscape through reiteration generating multi-multiple alignments (if that makes sense) using different gap penalties and considering the error between the alignments (highlighted below). The search strategy is key to the success of the output. Muscle 5 describes 'introducing variations (perturbations) into HMM parameters', I assume this is singularly with reference to super5. The parameter landscape might be quite small because it doesn't seem to be a directional search strategy in any way (there's no memory in it).
The optimisation algorithm, i.e. what is being observed is a metric called the 'pair-wise error (PWE) and pair-wise dispersion (PWD)', so what appears to be happening is the algorithm is attempting to minimise the number of gaps in an alignment and this would be compatible with your results.
Speed Super5 is breaking the alignment up into smaller alignment and solving each 'centroid' before bringing together the 'centroids'. This is an obvious strategy particularly in the era of multiple processors because it overcomes RAM bottlenecks, we are all doing that in any algorithmic context.
Overall The algorithm was not absolutely clear where PWE and the feedback loop is being deployed whether it was in the generation of centroids or bringing together the centroids and its relationship to the guide tree. There was a large reliance on BLOSSUM62 - which didn't quite make sense in the given context because the background was clear about the different substitution matrices and the impact on gap penalties.
Summary Static gap penalty parameterisation is antithesis of what super5 is about, that is why you cannot use them. The speed is clear its using parallel processing and subalignments - there are no surprises there, thats a common strategy and easy to code these days.
Final points

I am clear the ms has not been accepted and is available here https://www.biorxiv.org/content/10.1101/2021.06.20.449169v1.full#ref-10 . I don't want to comment further on the manuscript for several reasons, which I prefer to avoid explaining.
It is better to wait until the paper is published because it will likely be different.
The algorithm itself will unquestionably be subject to simulation studies by other authors.
Simulation studies were absent in the BioRXIV manuscript.
Muscle historically has been at the forefront of methods to increase sensitivity. If you read the Clustal Omega publication they describe the implementation of historical Muscle sensitivity improvements, in other words they must have been approved over the course of time.

Other authors are likely to investigate whether simulation conditions will result in super5 causing a 'false peak' scenario. That is the Achilles heal in any HMM.
There are certainly situations where super5 is useful in automation, that is without question. In some pipelines there is no option of checking the alignments within the pipeline because there are just to many and it stops the pipeline (human eyes the opposite of a data pipeline).
In practice optical annotation is used to overcome limitations post alignment, i.e. to resolve 'errors' in gap penalties and this is desirable because it prevents a phenomena of 'circular mathematics' which is a caution in any alignment method. The author does discuss this phenomena without calling this directly 'circular mathematics'. Any alignment algorithm can erroneously impose itself on the alignment and the metric of 'error' is extremely important in this regard. Optical analysis breaks circularisation. There is a very, very obvious solution to all this without using HMM at one extreme or 'optical analysis' at the other extreme, but one day at a time.
Example of HMM 'Circularisation' is the Achilles heal of my explanations, because I explain stuff, then reexplain it which more likely increases confusion and intersperse everything with waffle. However, I do follow a HMM style algorithm because next time I explain this it will be more honed and punchier. The hidden parameter being 'clarity' and the observed phenomena being the frequency of blank expressions, which is sought to minimise on each iteration. I do hope my search of the parameter landscape is a bit more directional than that used in Super5.

parameters looking at the documentation
-perturb SEED

It is important to do the same alignment using different random number streams. A given seed should result in the same gap pattern. If there are different gap patterns between different 'SEEDS', there's something suboptimal in the HMM and that is not good at all. If you think the Muscle 5 is implicit perfection, you don't need to check it - but its a trust thing. Its good the author made it explicit, i.e. you can set it.
The other parameters are:
-perm PERM
-consiters N
-refineiters N

The second most important parameter in my view is -perm - this is permutations on the UPGMA guide tree.  If different -perm values give different alignments thats not good. It must work reasonably well because the author will have checked it.

Thanks for the detailed explanation!


text:
Is SBOL-VOWS operational?
Synthetic Biology Open Language's Visual Ontology Webservice docs have example URLs but they all return a github pages 404.
EDIT:
The SBOL-VOWS has site has been fixed but the section on the SBOL website still has links with a 404.

would you mind adding the link in the docs that returns 404 github page? And perhaps mentioning what did you expect to see there instead? It just so people will be more easily understand/find your question if then have the same problem.

@KamilSJaron ah my bad. I just realized that SBOL website itself still has broken links!

Many thanks for raising this issue. The web service is now working. Only the glyph related web service calls were affected due to restructuring folders at the SBOL-Visual GitHub repository. I now revived the glyphs using the sbol-visual-ontology GitHub repository. Over the summer, we will develop the web service further so that you can request different versions of glyphs. Currently, the web service works with the SBOL Visual v2.3 glyphs and we will provide support for SBOL Visual v3 glyphs soon.
The web service is available at:
https://vows.sbolstandard.org/
Please let me know if there is anything I can clarify.

Hello - I've just noticed the links on the SBOL website are still broken.


text:
Combining dataTracks in Gviz by averaging or subtraction
I would like to average (and maybe also additional operations such as subtraction) 2 or more DataTracks and plot the result using the Gviz (https://bioconductor.org/packages/release/bioc/html/Gviz.html) package.
Say I have:
chr="chr2"
gen="mm10"
from=1000
to=5000

dTrep1 <- DataTrack(range = path/to/rep1.bw, genome = gen,
             name = "bigWig",chromosome = chr)
dTrep2 <- DataTrack(range = path/to/rep2.bw, genome = gen,
             name = "bigWig",chromosome = chr)

Instead of:
plotTracks(list(dTrep1,dTrep2),type="polygon",from=from,to=to)

I would like to:
plotTracks(mean(dTrep1,dTrep2),type="polygon",from=from,to=to)

Or:
plotTracks(dTrep1-dTrep2,type="polygon",from=from,to=to)

See posted also at: https://support.bioconductor.org/p/9137211/
Thanks.

It might help to post this on the BioConductor support forum and add a link to that post here, either as a comment or by editing your post and adding it in the body.


text:
How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
I would like to download a number of SARS-CoV-2 fastas from NCBI/Genbank using a web endpoint. I looked on their website but it's confusing.
I have the accession numbers, but can't figure out how to use the NCBI/Genbank Entrez API.
This is how I would do it using cov-spectrum: https://lapis.cov-spectrum.org/open/v1/sample/fasta?genbankAccession=MZ039713,MZ039682,OV522182
Is there something similar for NCBI/Genbank?
There must be a way, I just can't figure it out. The whole entrez thing is overwhelming to me, all I want is a simple API, like covSpectrum offers, for example. I scanned this but didn't help: https://www.ncbi.nlm.nih.gov/sars-cov-2/
E.g. how could I download all these SARS-CoV-2 nucleotide sequences with accession numbers: MZ039713,MZ039682 and OV522182 using a URL?
Some command line tools are advertised, but I don't want to have to use a CLI tool. Just a simple URL, does that not exist? CLI tools: https://www.ncbi.nlm.nih.gov/datasets/docs/v1/quickstarts/command-line-tools/
I also don't want to do any clicking through a web interface, just a URL so that other web apps can download it behind the scenes.

Could you please let us know which ones?

@M__ what do you mean by "which ones"? For example fastas for these accession IDs: MZ039713, MZ039682 and OV522182 like displayed here: https://www.ncbi.nlm.nih.gov/nuccore/OV522182

Thanks I'll put a response together. Can't do it immediately, but will do today

You can get the data from the API of the European Nucleotide Archive (ENA) which is synched with NCBI/GenBank:
https://www.ebi.ac.uk/ena/browser/api/fasta/MZ039713,MZ039682,OV522182

You can try the NCBI Datasets API. We have an endpoint that would allow you to download virus sequences by accession  The only thing to be aware (based on your list of accessions) is that our API needs the accession version to return results. In the case of your examples, they would need to be MZ039713.1, MZ039682.1 and OV522182.1.
I hope this helps!

Thanks! Can you provide an example query? The API documentation is a bit overwhelming to me. Also, am I correct in assuming you need a private API key?

You don't need a private API key.
Using the accessions you provided as an example, you can create a json file (accessions.json):
{"accessions":["MZ039713.1", "MZ039682.1", "OV522182.1"]}
Then, you can use the following command to download everything using the API:
cat accessions.json | curl -H "Content-Type: application/json" -X POST --data-binary @- "https://api.ncbi.nlm.nih.gov/datasets/v1/virus/genome/download" --output file.zip
Let me know if you run into any issues.

In the efetch.fcgi URL, you can set rettype=fasta to download the nucleotide sequence or rettype=fasta_cds_aa to download protein coding sequences:
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta&id=NC_045512.2'|head -n3
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta_cds_aa&id=NC_045512.2'|head -n3
>lcl|NC_045512.2_prot_YP_009724389.1_1 [gene=ORF1ab] [locus_tag=GU280_gp01] [db_xref=GeneID:43740578] [protein=ORF1ab polyprotein] [exception=ribosomal slippage] [protein_id=YP_009724389.1] [location=join(266..13468,13468..21555)] [gbkey=CDS]
MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLVEVEKGVLPQLEQPYVF
IKRSDARTAPHGHVMVELVAELEGIQYGRSGETLGVLVPHVGEIPVAYRKVLLRKNGNKGAGGHSYGADL

There are three ways to do this,

If you click on this link Your NCBI examples. Its done.
You would then just click "download"
Click "next", click "next"

The sequence then downloads.
Essentially all I did was in the website click on the left hand navigation bar and clicked "accession number" and type in MZ039713 MZ039682 OV522182.

The second way is to go to https://ncbi.nlm.nih.gov/nuccore
Type type in MZ039713 MZ039682 OV522182
Click on the "Summary" and click on "Fasta"
Click on "send to" and click on "File"

I'm not certain is point 3 immediately above is required.
The final approach is to use efetch to automate the download, but that requires downloading the NCBI's efetch suite, which is then flipped through on a bash loop.
efetch -db nuccore -id MZ039682 -format fasta

The link for looking to automate is:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide&VirusLineage_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20
Just to highlight
MZ039713%20MZ039682%20OV522182%20
The API is here, but I don't use them.

Thanks! I clarified that I would like the download link to be programmatic = API, not needing to do any clicking in a web app.
So do I understand it correctly that there really is no simple API to download these fastas? I need to download some efetch suite to do programmatic = non-interactive downloads?


text:
Correspondence of SARS-CoV-2 annotations (Nextstrain clades - Pango lineages)
Is there any annotation file providing correspondence between Nextclade and Pangolin variant nomenclature/annotations, to annotate some SARS-CoV-2 genomes from Gisaid with both these?
For the moment I could only find:

Nextclade's clade definitions
Pangolin's lineage designations

What I'm looking for is an annotation as follows:
lineage    clade
AY.43      21J
AY.4       21J
...        ...

The idea is to generate a table as follows:
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

The way I proceed now is the following:

Parse the Pangolin's lineage designation file to find a target genome
Look for the target genome in Gisaid (for quality check)
Parse the metadata associated to the selected genome to get the lineage for verification (only Pangolin lineage is reported)

These steps allow me to build the first part of the table:
genome                                 lineage
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43
...                                    ...

Then:

Paste the genome to Nextclade, wait for the analysis to be completed
Report the calculated clade

With this information, I can then build the table:
genome                                 lineage    clade
hCoV-19/Germany/BW-RKI-I-195742/2021   AY.43      21J
...                                    ...        ...

The idea behind this is that both lineages and clades can be used for "tagging" a genome, depending on the aim (a given clade e.g. 21J includes many lineages e.g. AY.43 and AY.4 in the example above) and that I don't enjoy manual and error-prone procedures;) Thanks for your input.

The GISAID nomenclature via Nextstrain and overall tree is presented here
(Crossposted on Nextstrain's discord forum)

The problem is that there's no one to one mapping between pango lineages and Nextstrain clades. One could generate the mapping from the Nextclade reference tree, there most tips have a pango lineage and a Nextstrain clades, and almost all pangos lineages appear in that tree. Other option is to use the open data metadata TSV from Nextstrain, same info.

A comment from a Nextclade's user satisfies my actual needs, I will post it here for future reference:

Hi - one possible source is the metadata file for UCSC’s daily-updated tree of public SARS-CoV-2 sequences:

https://hgdownload.soe.ucsc.edu/goldenPath/wuhCor1/UShER_SARS-CoV-2/public-latest.metadata.tsv.gz

It’s a big file – almost 3 million lines – because there are about 3M sequences in GenBank & COG-UK at this point and that number will probably grow. The last two columns of the file (10th & 11th) are the Nextstrain clade and Pango lineage assigned to each sequence based on its placement in the UCSC/UShER tree and the tree’s annotation of where each clade/lineage starts.

This command will make a file relating each Nextstrain clade to at least one lineage (often to a lot of lineages, especially for Delta):

gunzip -c public-latest.metadata.tsv.gz | cut -f 10,11 | sort -u > cladeToLineages.tsv

Some Pango lineages have no public sequences (only GISAID sequences, which can’t be shared at a public URL), so those lineages will be missing from that file. Some older lineages are not monophyletic on the UCSC/UShER tree, so those lineages are not annotated and will be missing. (In total, ~200 / ~1500 are missing.)

This information is now easily available in JSON format here: https://github.com/corneliusroemer/pango-sequences/blob/main/data/pango-consensus-sequences_summary.json. It's one big dictionary of dictionaries with Pango lineage names as outer keys and corresponding Nextstrain clades accessible via inner key:
{
 "A": {
  "lineage": "A",
  "nextstrainClade": "19B",
  "nucSubstitutions": [
   "C8782T",
   "T28144C"
  ],
  ...
 },
 "A.1": {
  "lineage": "A.1",
  "unaliased": "A.1",
  "nextstrainClade": "19B",
  ...
 },
 ...
}

It's straightforward to convert this to a TSV mapping using just minimal scripting.
Beware that not all Pango lineages pre-Omicron are contained therein. Small ones are missing currently.
Disclaimer: I created the repo.

Remember not to omit the 'clarification' that in addition to Pango lineages, GISAID clades, Nextstrain clades there is also the WHO nomenclature. We could designate the Bioinformatics StackExchange monophyly, surely there is room for one more nomenclature?
Here's the WHO's response copied and uploaded from here

Just  to mention NCBI has adopted the Pango lineage nomenclature available on all of its sequences available via this link
I thought that Omicron was Pango lineage B12, or perhaps B.1.2? Oh I see it has been reclassified following further sequence data.


text:
Download Addgene database
Is anyone aware of a way to download all vector sequences from addgene at once? As a non-profit, I'd assume they shouldn't be opposed to making this data available, but it doesn't seem to be possible. Might anyone have scraped it at some point?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

That is a very good question for Addgene. You could contact them directly.


text:
Smith-Waterman traceback with large gaps
I am currently trying to understand how the traceback algorithm is supposed to work for the smith-waterman algorithm as my current understanding breaks down in case of a large alignment gap.
Assume the sequences
3GA12CA
GA56CA

The obvious best alignment would be (or some other variation thereof)
3GA12-CA
-GA--6CA

However, assuming an affine gap penalty function of x => x + 2 and a binary similarity score (i.e., if equal then 1 else -1) would lead to a scoring matrix
    3 G A 1 2 C A
  0 0 0 0 0 0 0 0
G 0 0 1 0 0 0 0 1
A 0 0 0 2 0 0 0 0
6 0 0 0 0 0 0 0 0
C 0 0 0 0 0 0 1 0
A 0 0 0 1 0 0 0 2

When doing the traceback I would start with the largest value towards the lower right corner of the matrix

(6,6)
(5,5)
Now I do not know how to find the next best alignment (in this example at (3,2)) as neither going up nor going to the left nor going diagonal is going to lead to a score greater 0 meaning I have no next best value to go on.

How is the traceback function supposed to continue here?
Additionally, I read that the traceback should start at the highest score in the matrix itself. However, if we would extend the first subsequence match length by one, this would mean we miss the right alignment option completely.

Before we get started, remember that mismatches are often preferred over gaps, so the "obvious" best alignment would treat one of the 12 in seq1 as a mismatch and the other as a gap, not introduce gaps in both sequences unless your mismatch penalty >> the gap open penalty (which makes no sense)

I'm still trying to understand your problem statement. Affine gap penalties make expanding gaps less expensive than opening new ones, they do not make gaps cheaper than mismatches. "Gaps" are biological sequence insertions/deletions, and those are serious changes, you do understand that, right?

This question is based on a misunderstanding. The Smith-Waterman algorithm calculates the best local segment. You have to start at the $max(H_{i,j})$ point and run the traceback as described. It will (most likely) return only a local segment that, however, has the maximum value in terms of alignment score.
To get a global alignment, another algorithm has to be used. Needleman-Wunsch would be an example.
The difference between local and global is that a local alignment only alignes a subsequence of the provided sequences whereas a global alignment gives the best alignment of the whole sequences.


text:
Remove from Multi-FASTA by Sequence ID
I want to remove sequence VRE32514 – it doesn’t belong and thus is the reason it lacks additional metadata. However I tried implementing this code from a similar question:
awk 'BEGIN{while((getline<">VRE32514")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' *.fasta | grep -v ">VRE32514" *.fasta

With no success and the proposed solution in R using seqinr doesn't seem to work either for myself.
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("test.fasta")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), "|"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c(">VRE32514")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fasta")

Below is an example of what my Multi-FASTA appears as and the desired output:
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT
>VRE32514
AATCAGATTGTTCGCCATGAGCTGCTGAGAATGATATCATGGAAGGGCGGCATCGAAACA
GGCTTTAAATTAAGTGTAGGCAAGAACTATAAGTTTATTGAAAGGTATATATCCGAGGAT

Desired Output:
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

Here's a solution using Python:
from Bio import SeqIO

sequences = SeqIO.parse('test.fasta', 'fasta')
filtered = [seq for seq in sequences if 'VRE32514' not in seq.id]

with open('output.fasta', 'wt') as output:
    SeqIO.write(filtered, output, 'fasta')
```


text:
Methods to work directly with bigbed file?
I am now trying to find intersect genes of two files: one is .bigbed file (very large, 117GB), and the other is a .bed file.
At first, I tried to convert bigbed file to bed file and then applied bedtools. I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server. Thus I planned to convert the .bb file on my own computer and then upload the converted .bed file to the server, and then apply bedtools. Firstly I tried to convert some smaller .bb data sets (< 10 GB) to .bed and it worked slowly but well on my computer. However, when I tried the 117GB data there was an error:
Invalid argument, Error reading 9188986873 bytes.
My code is :
bigBedToBed JASPAR2022_hg19.bb hg1 9.bed
Thank you in advance!

Now the bug is solved, and I got many other good suggestions. Below you can find many methods to work with the bigbed file, directly working with it or converting it to bed file.

You asked

Does anyone know any tools that can intersect two bigbed files

...but it appears you really wanted to subset a single very large bigbed file based on the contents of a much smaller bed file.
I disagree with your summation "there is not a very efficient way to work directly with the bigbed file".  I see two options.
The first is to repeatedly call bigBedTobed once for each row in your .bed file.  If you have multiple cores, this can be parallelized using GNU parallel.  The following does that, emitting each overlapping result with additional column being the index into your .bed file from which the result stems.
cat your.bed | parallel --lb -j --colsep '\t' bigBedToBed -chrom={1} -start={2} -end={3}  your.bb  /dev/stdout '|' sed -e 's/$/\\t${PARALLEL_SEQ}/'

Depending upon your application, there is a second approach which uses R/BioConductor's rtracklayer::import.bb which can be parallelized across (chunks of) your .bed file.
I considered these approaches in what I expect is an almost identical situation:
My bigbed holds genome-wide transcription factor binding site predictions from JASPAR Genome Browser tracks as produced by wassermanlab / JASPAR-UCSC-tracks.
My bed file holds the loci of ~23K putative enhancers as determined by an analysis of regions of accessible chromatin from ATAC-Seq and various ChIP marks.
In my case I wished to produce a sparse Matrix with one row for each putative enhancer, and one column for each JASPAR motif (which will form the basis of downstream enrichment analyses).
bbf<-BigBedFile('path/to/your.bb')
gr<-import.bed('path/to/your.bed')
## I found I had to harmonize seqlevels (e.g. chromosome identifiers) to make this work.  YMMV
gr<-renameSeqlevels(gr,paste0('chr',gsub('^(K.*)\\.(\\d+)$','Un_\\1v\\2',seqlevels(gr))))
mm<-motifMatrix(gr,bbf)

which results in a sparse Matrix, viz:
> rownames(mm)<-as.character(gr) # optional.
> mm[1:5,1:5]
5 x 5 sparse Matrix of class "dgCMatrix"
                 MA0002.2 MA0003.4 MA0004.1 MA0006.1 MA0007.3
chr1:11156-12180       14        .        2       28        2
chr1:17836-20023       25        .       27       63        2
chr1:22495-24367       22        1       22      137        2
chr1:27533-28038        2        2        4        7        .
chr1:36549-37079        6        .        2       41        .

This depends upon the following definition:
motifMatrix<-function(gr,bbf,motifID=c('name','TFName')[[1]],mc.jobs=100) {
  ## produce a sparse Matrix with one row for each region in
  ## GenomicRanges gr, and one column for each JASPAR motif indicated
  ## as overlapping <gr> accorinding to BigBedFile <bbf>.  Default to
  ## using the 'name' column as the motif identifier, but allow for
  ## possibly using another, such as 'TFName' (as used by 2022 version
  ## of
  ## https://jaspar.genereg.net/genome-tracks/#ucsc_tracks). Parallelize
  ## into mc.jobs queries, defaulting to 100, which was found suitable
  ## for ~25000 loci in <gr> against JASPAR matches to zebrafish
  ## danRer11 genome running on a few dozen cores.  In general,
  ## knowing how "big" each parallel job *should* be depends on
  ## available RAM and cores, and number of loci in <gr>.
  ##
  ## AUTHOR: malcolm_cook@stowers.org
  ##
  ## NB: you must install a version of `rtracklayer` which addresses
  ## issue [import\.bb returning incorrect ranges & results as
  ## compared with kent tools
  ## bigBedToBed](https://github.com/lawremi/rtracklayer/issues/59)
  ## (such as by `install_github('lawremi/rtracklayer')`)
  si<-parallel::splitIndices(length(gr),mc.jobs)
  mm<-mclapply(si,function(i) {
    gri<-gr[i]
    m<-import.bb(bbf,selection=BigBedSelection(gri, colnames = c(motifID)))
    fo<-findOverlaps(gri,m)
    xt<-
      xtabs(~.,data=list(locus=i[queryHits(fo)],motif=mcols(m[subjectHits(fo)])[[motifID]]),sparse=TRUE)
    xt
  })
  mm<-
    ## rBind.fill(mm) ## NOT!  It does NOT preserve dimnames and
    ## returns matrix when given Matrix.  Issue reported:
    ## [[https://github.com/cvarrichio/Matrix.utils/issues/5][rBind.fill
    ## on a list of Matrix arguably should return Matrix but returns
    ## matrix #5]]
    Reduce(rBind.fill,mm) # which does preserve dimnames and return Matrix
  mm
}

NB: If you choose this approach, you MUST install a version of rtracklayer which addresses issue the import.bb returning incorrect ranges & results as compared with kent tools bigBedToBed (such as by install_github('lawremi/rtracklayer'))

Yes Malcook, this is exactly what I want at first! Thank you so much!

glad to help - would you consider to change the accepted solution to mine?

This is the accepted answer for sure @BinhuanSun

My 2p, a little late and maybe not too useful:

I didn't run the bigbedtobed program on the server because I am not allowed to install the bigbedtobed program on our server.

Presumably, you cannot install in /usr/local/bin or some other root directory because you don't have admin rights. This is the common setup on shared servers. However, many programs can be installed in the user's home directory without admin rights. Even better, install conda and mamba first (no need to be admin) and use that to install other programs.
The problem with docker is that some system administrators are concerned about security and refuse to install it (and docker does need admin rights). Whether such concerns are founded, I don't know...

I tried to convert bigbed file to bed file and then applied bedtools, but I found it was very slow and couldn't make it.

bedtools intersect supports the -sorted options that should speed things up quite a bit - Have you applied it? Your huge bigbed files should be already sorted after you convert to bed.

Hi Dariober, yes later I found I could install programs in my home directory. Sorry for this super easy problem because I am not familiar with the Linux server. I just checked my converted bed file, and it is indeed sorted, so I think it is just because the converted bed file is too big (560GB).

No matter how you approach this problem 110+ GB is big:
Your options are:

Break up the file and perform the intersection on a chunk by chunk basis
Ask your system administrator to install the program
Load it into Docker, your systems administrator will understand this. Installing docker is a single line of code (which I forget). Docker is a container which protects the system from the program and no systems administrator in this day and age would refuse a request for Docker installation.

Your.systrm administrator might be concerned about the file sizes involved, but I think a better dialogue with them is the right way forward here.

the first option is maybe a good one, you can select individual chromosomes at a time with bigBedToBed

Thank you for your suggestions! I will try the first one on my computer first :)

@ColinD Yes, I will try. Thank you :)

Hello M, because today is Sunday so I first tried the first method you suggested. I tried to break up the file using argument -Chrom of the bigBedtobed program: bigBedToBed -chrom=chr1 JASPAR2022_hg19.bb test.bed. However, I still got the same errorInvalid argument, Error reading 9188986873 bytes. I then tried a relatively smaller data (66GB) but got the same error. It seems that on my computer the bigBedtobed program only works on data that is smaller than 10 GB. I still have over 600 GB of drive memory, so I think this may due to the lack of RAM?

If its been written in C (I honestly dunno) a core-memory dump is a common error. However, if you > 10 chromosomes each file should be smaller than 10 GB? If its Drosophila then you get the same error. Could you check the file size for each intersect?

Hello M, Thank you for your quick reply! Actually, it is a Genomic track and TFBS prediction for human genome, and it is a whole file. I think after breaking up every file will < 10GB, but the thing is I don't know how to split the .bb file directly and then convert them to the .bed file :(. The website for the file I am working with is here: link. The file name is JASPAR2022_hg19.bb.

I think this should be a separate question. It wouldn't be hard to modify the file to allow breakup the file given the chromosome tag is a recognised tag for this purpose as @ColinD has stated. Perl pie will modify massive files with ease (if it was simplistic).

Yes M. I think you are right, and today I asked my supervisor to help me apply for an account on another server that installed those UCSC programs, and I hope it can work on the server. Thank you again for the patient reply :)

The question I asked a couple of days ago has been solved. Here I want to make a summary of this and hope it can help people who also encounter such problems later.
First, about the way to process the bigbed file. I didn't find a very good way to work directly with bigbed file at that time, so I still firstly converted the bigbed file to the bed file using bigBedTobed program and then apply bedtools. The bigBedTobed program can be downloaded through this link: UCSC programs, and it also contains many useful programs from  UCSC genome browser. I finally converted my data (117 GB bigbed file) on the server, and it took ~2.5h to convert it to bed file. After conversion, the converted bed file is about 550 GB. And it took ~4h to intersect with a 4 Mb bed file using bedtools.
Second, the bug that I encountered when I tried to convert the bigbed file on my own laptop. The code I used is bigBedToBed JASPAR2022_hg19.bb hg1 9.bed and it has error Invalid argument, Error reading 9188986873 bytes. I reported this error to the UCSC genome browser group and they said it is a MacOS specific bug, reaching a default data limit. They have solved it in their development environment but it still needs some time to reach the public site. The current solution for this is to use a URL to the file instead of a local file allows the program to read appropriately sized pieces of data at a time(-udcDir= in bigBedToBed).
Finally, this is my first time asking questions here and I really appreciate those good suggestions from @M__ and @ColinD.

Renew: Methods from @malcook can work directly with the bigbed file! For a very large bigbed file, it is indeed a more efficient way.

Hi nice to hear from you @BinhuanSun. Cool glad its worked out and thanks for the feedback.

MacOSX bug .. thanks thats interesting, I encountered something similar. They are nice machines, but it ain't a cool way to do things IMO.


text:
Time complexity of the Neighbors function
I implemented the following find_neighbors_with_expected_hamming_distance function that generates all k-mers of Hamming distance at most d from the given Pattern.
from .Hamming_Distance import find_hamming_distance

dp_map_that_stores_all_neighbors_with_expected_hamming_distance = {}

unique_nucleotides = {'A', 'C', 'G', 'T'}

def suffix(pattern):
    if len(pattern) <= 1:
        return None
    else:
        return pattern[1: (len(pattern))]

def first_symbol(pattern):
    if len(pattern) == 0:
        return None
    else:
        return pattern[0]

def find_neighbors_with_expected_hamming_distance(pattern, expected_hamming_distance):
    if expected_hamming_distance == 0:
        return {pattern}

    if len(pattern) == 1:
        return {'A', 'C', 'G', 'T'}

    if (pattern, expected_hamming_distance) in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        return dp_map_that_stores_all_neighbors_with_expected_hamming_distance[(pattern, expected_hamming_distance)]

    neighborhood_set = set()

    suffix_neighbors_set = find_neighbors_with_expected_hamming_distance(suffix(pattern), expected_hamming_distance)

    for suffix_neighbor in suffix_neighbors_set:

        if find_hamming_distance(suffix(pattern), suffix_neighbor) < expected_hamming_distance:

            for nucleotide in unique_nucleotides:
                neighborhood_set.add(nucleotide + suffix_neighbor)

        else:
            neighborhood_set.add(first_symbol(pattern) + suffix_neighbor)

    if (pattern, expected_hamming_distance) not in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        dp_map_that_stores_all_neighbors_with_expected_hamming_distance[
            (pattern, expected_hamming_distance)] = neighborhood_set

    return neighborhood_set

I want your suggestion on my time complexity calculation:
For a given pattern of length N and D as maximum expected hamming distance, I calculated time complexity as -
$\sum_{n=1}^N$ $\sum_{d=1}^D$ $n\choose d$ * $3^i$
Reasoning: For any length n, we will choose any d letters (nucleotide in this case) to be replaced with the other 3 nucleotides that is different from the nucleotide present at that particular position.
For example: Pattern = ATGCAT and D = 2, then suppose for one of the case, I select 2nd and 4th positions to be replaced.
Then, at the second position (current nucleotide = 'T'), I can substitute any nucleotide from the set {'A', 'C', 'G'} and
for fourth position (current nucleotide = 'C'), I can substitute any nucleotide from the set {'A', 'T', 'G'}.
Kindly tell me if this approach is correct?
EDIT:
I have provided the pseudocode for the Neighbors function as follows:
Neighbors(Pattern, d)
    if d = 0
        return {Pattern}
    if |Pattern| = 1
        return {A, C, G, T}
    Neighborhood ← an empty set
    SuffixNeighbors ← Neighbors(Suffix(Pattern), d)
    for each string Text from SuffixNeighbors
        if HammingDistance(Suffix(Pattern), Text) < d
            for each nucleotide x
                add x • Text to Neighborhood
        else
            add FirstSymbol(Pattern) • Text to Neighborhood
    return Neighborhood

The symbol • means concatenation of 2 strings. For example, 'S' • 'V' = 'SV'.
If we remove the first symbol of Pattern (denoted FirstSymbol(Pattern)), then we will obtain a
(k − 1)-mer that we denote by Suffix(Pattern).
Illustration:
Given pattern = 'ACG'
Expected Hamming distance (D) = 1
Answer = {'CCG', 'TCG', 'GCG', 'AAG', 'ATG', 'AGG', 'ACA', 'ACC', 'ACT', 'ACG'}
Note that the Answer is a set of all strings who have hamming distance of 1 with respect to the given pattern 'ACG'.

I'm confused by the "expectation" part. Isn't this deterministic?

Also, I'd recommend adding a pseudocode explanation of your approach in addition to the python code.

@Throckmorton, thank you very much for having a look. I have added the pseudocode with an illustration.

Let's first establish a lower bound on any solution:
Given a hamming distance $D$ and an input pattern $P$ of length $N$, there are $\sum_{k=0}^D{N \choose k} 3^k$ other strings within hamming distance $D$ of $P$. Of course in the worst case, when $D=N$, there are $4^n$ possible strings.
We can also get a closed bound on the number of solutions. We can show that the sum of the first $D$ binomial coefficients for a fixed $N$ is bounded by ${N \choose D} {N-(D-1) \over N-(2D-1)}$ as shown here. Therefore we have the following closed upper bound on the number of strings within hamming distance $D$ of some pattern of length $N$
$$
\sum_{k=0}^D{N \choose k} 3^k \leq {N \choose D}  {3^D (N-(D-1)) \over N-(2D-1)} \leq \left(\frac{3N}{D}\right)^D {N-(D-1) \over N-(2D-1)}
$$
However, for your approach, you are going through recursively, so when the subproblem is on strings within $D$ of the suffix of length $n$, you have $\sum_{k=0}^D{n \choose k} 3^k$  of these. Given that you have $N$ subproblems (suffixes), at first glance it does indeed look like you have roughly $\sum_{n=1}^N\sum_{k=0}^D{n \choose k} 3^k$ iterations of the innermost for loop. However, in that for loop, you are also computing the Hamming Distance function, which is $O(n)$. This means that you have roughly
$$
\sum_{n=1}^Nn\sum_{k=0}^D{n \choose k} 3^k
$$
computations.
I would suggest finding a way to get rid of having to compute the hamming function. Since you are constructing the sequences yourself, you can keep track of each sequences hamming distance from the pattern during construction.

Excellent! I missed the complexity for computing Hamming distance. Also, you can correct the first sentence - "Given a hamming distance  D"

Thank you very much @Throckmorton.