title: Affinity maps and forcefields

text:
# Question
Title: Affinity maps and forcefields
I know that before performing docking with AutoDock Vina, we have to compute "vina maps" like below

```
v.set_ligand_from_file('1iep_ligand.pdbqt')
v.compute_vina_maps(center=[15.190, 53.903, 16.917], box_size=[20, 20, 20])

```

but I do not quite understand what those maps exactly are. I found from the manual that

> ...and then compute the affinity maps for each ligand atom types accroding to the Vina forcefield.

but again, what are affinity maps and Vina forcefield. I checked the source code to see that it creates some kind of a grid which will be then used for scoring. What kind of information does this grid data type contains? I would greatly appreciate any help.

# Answer
Docking is a computational calculation in the field of **molecular mechanisms**. It can be argued that whereas it is true one can 'wing it' without knowing the basics, it is best to know them —akin to the fact that it is possible to drive a car fine without knowing the road code or how a car works: many do for sure but it is not a good idea as something will go wrong.

So I would suggest having a gander at the very least at:

So a forcefield is a sum of different terms that mathematically describes the forces acting between atoms. A thing to note is that these terms may be have some basis in chemistry or some may be just fudge factors to make the data fit.

Regarding the atom types, in a full atom model, each atom is represented as a particle (or node in graph theory speak) and the bonds are fixed (edges). An element behaves different depending on its hybridisation and charge etc. Therefore, each atom is classified not as an element, but as an atom type. Say, nitrogen in a ring and amine are generally different atom types. With parameterising a small molecule an atomtype may not be a perfect fit.
> 1 votes

---
Tags: docking
---


# Question
Title: Is this a complete bacterial genome? (CheckM)
I've been working with a bacterial genome assembly. The initial stats look pretty good (3 main pieces, depth ~50x \[20x - 100x\] if I calculated it right). I'd expect this to be a virtually complete genome apart from some low-complexity regions.

However, when I use CheckM to estimate completeness / redundancy, I get ~99% completion, \<2% redundancy. (out of 579 expected single-copy marker genes, 563 were present once, and 11 were present twice). I am especially concerned about the missing genes.

I can think of several possibilities: a) Limits of CheckM. This software uses a sample of 82 genomes to select marker genes for this phylum: however, there are no available genomes in the same genus as the one I'm working with. It might be possible that this is natural variation.

b) The "single-copy" genes present twice might be due to misassembly. But while the sample was not axenic, the organism of interest was majoritary, and long reads were used as well, so this shouldn't be an issue?

c) The genes were somehow missed by sequencing. (but the depth is pretty good?)

I don't know how likely a) is, so it's difficult for me to assess whether what I have is a normal result, or something that requires additional sequencing in order to make sense. I would appreciate your thoughts!

# Answer
> 2 votes
See this table from the CheckM table to put this number in context. In general, your numbers sound excellent, even compared to isolate genomes.

I'd recommend for context taking some finished microbial genomes (E. coli K12 etc) and running CheckM. I bet you'd be surprised what you find.

I think that you can probably not worry about (c) as a possibility. If you are very worried about (b), you could map reads against your assembly, filter out the ones that don't map well, and assemble those to see what you find. Long reads suggest that your genome is in good shape. You can of course do manual curation of your genome using tools like Anvi'o, or by scaffolding/gap-filling your assembly with long reads or an orthogonal strategy like optical mapping or Hi-C.

I think that the influence of (a) is significant. Overall it sounds like you have a reasonable genome.

# Answer
> 2 votes
Don't read into the CheckM result too much. Genes having one copy in some species may have multiple copies in other species and genes present in some species may be absent in others. This may even happen between strains, let alone between species or between genera. It is rare to see 100% checkM completeness and 0% contamination. A common standard for "near complete" assembly is \>90% completeness and \<5% contamination. Your assembly is far above this standard.

Meanwhile, beware that checkM only evaluates single-copy genes but multi-copy genes are usually harder to assemble. There might be complex multi-copy genes in the three assembly gaps. Without getting a circular contig, it is hard to say for sure your assembly is really complete.

---
Tags: genome, quality-control, genome-sequencing
---


# Question
Title: Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
I implemented Hirschberg's algorithm in python and used the wiki example to verify correct implementation given the scoring parameters and sequences:

```
Let X = AGTACGCA
    Y = TATGC
Del(X) = -2
Ins(Y) = -2
Sub(x,y) = -1 for mismatch +2 for match

This aligns to:

AGTACGCA
--TATGC-

```

Which matches the wiki output. I also printed the recursion and it matched the wiki recursive partitioning too. Then I tried it on a harder problem and got a terrible result. It was a mitochondrial Neanderthal DNA sample of about 500 nucleotides and a complete human DNA mitochondrial sample of 16504 nucleotides. The output of trying it on BlastN is here:

https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr\_Query\_40035

My results look something like this where the top is neanderthal and the bottom is the full human genome. It appears though it just matches each character as soon as it finds on.

```
---C-CA-----A--------------------G-----T----AT---T--G--A----C-TC--------A--C-C-C-AT--CA-----A--C-----A-----A-C-C-----G-C-CA-T----GT-----ATT--T-C------------------G----TAC----A-T-T-AC---------T----------G------------------C--------CA--------------G-----C--CAC--C-A-T---------GA-AT-AT---------T-----------------------G-T-------ACAG---T-A--C-CAT-------AA--------------------T----------TA-C----TT----------------G-----A--C--T---A-C---C----------T--------------------------------------------------------G----T---A---A---------------T---A-C---ATA-----AA--AACC------------------------T--A---A--T--CC-----A---CA-T-CA---A-----------C---C-C-C--C-CCCC-----C-------------CC-A---T--------GCT-T-A---------CA-A-GCAAGCA-C--------AG------CA------A--TCA--AC---C--------------T--------------T-CA----A---C--T--G--T--C-A---------------------T-A----C----A----T---C-AA-------CTA--C-A-ACT--CC--A---------AA----G----A-C-ACC-C--T-------T-A-C--A---C--------CC-----A---------------C------T----A-----G-------------GA-T---A------TC-A----ACA-AA------C--C-----T------A-C----C----CAC-C-----CT-----------TG--A-------C---A----G-T-A-C---A----T-A---G------CA-CA----T-------A-AA------G---------T-----C--A-------T-----T--T--A---CC-----G-------T-------A--C-ATA------G--CA-C---A-----T-T-----A---TA------G---TC---A--AA----T--------C--C--------C-----T-------T-----C-T-----C--G------CCC---------C----------C-A--T-----------G--------GAT-G-----A-----C--------C--------------C------C-------------C---------C-------T----C----A------G------A---T----A-G-----------GG---------------------G--T--------------C----C----C-T--T--G---A---------------------------------------------------------------------------------------------------------------
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAATGTCTGCACAGCCACTTTCCACACAGACATCATAACAAAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCAATACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCAAGCACGCAGCAATGCAGCTCAAAACGCTTAGCCTAGCCACACCCCCACGGGAAACAGCAGTGATTAACCTTTAGCAATAAACGAAAGTTTAACTAAGCTATACTAACCCCAGGGTTGGTCAATTTCGTGCCAGCCACCGCGGTCACACGATTAACCCAAGTCAATAGAAGCCGGCGTAAAGAGTGTTTTAGATCACCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGGCTTTAACATATCTGAACACACAATAGCTAAGACCCAAACTGGGATTAGATACCCCACTATGCTTAGCCCTAAACCTCAACAGTTAAATCAACAAAACTGCTCGCCAGAACACTACGAGCCACAGCTTAAAACTCAAAGGACCTGGCGGTGCTTCATATCCCTCTAGAGGAGCCTGTTCTGTAATCGATAAACCCCGATCAACCTCACCACCTCTTGCTCAGCCTATATACCGCCATCTTCAGCAAACCCTGATGAAGGCTACAAAGTAAGCGCAAGTACCCACGTAAAGACGTTAGGTCAAGGTGTAGCCCATGAGGTGGCAAGAAATGGGCTACATTTTCTACCCCAGAAAACTACGATAGCCCTTATGAAACTTAAGGGTCGAAGGTGGATTTAGCAGTAAACTAAGAGTAGAGTGCTTAGTTGAACAGGGCCCTGAAGCGCGTACACACCGCCCGTCACCCTCCTCAAGTATACTTCAAAGGACATTTAACTAAAACCCCTACGCATTTATATAGAGGAGACAAGTCGTAACATGGTAAGTGTACTGGAAAGTGCACTTGGACGAACCAGAGTGTAGCTTAACACAAAGCACCCAACTTACACTTAGGAGATTTCAACTTAACTTGACCGCTCTGAGCTAAACCTAGCCCCAAAC

```

The neanderthal sample matches at the end of the human mitochondrial dna. I suspect the gap score penalty is far too costly to insert enough gaps to reach the end of the alignment but I set blast to a linear gap penalty. How do I overcome this? Here is the implementation code for my program:

```
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

# NWScore
def NWScore(X, Y):
    # Initialize NM Scoring Matrix. Array needs to be the length of X +1 and Y + 1
    score = [[0 for i in range(len(Y)+1)] for j in range(len(X)+1)]
    for j in range(1,len(Y)+1):
        # Insertion scoring penalty column
        score[0][j] = score[0][j-1] + ins("a")
    for i in range(1,len(X)+1):

        # Delete Penalty Row
        score[i][0] = score[i-1][0] + dele("a")
    for i in range(1,len(X)+1):
        for j in range(1,len(Y)+1):
            # Scoring Matrix
            scoreSub = score[i-1][j - 1] + sub(X[i-1], Y[j-1])
            scoreDel = score[i-1][j] + dele(X[i-1])
            scoreIns = score[i][j - 1] + ins(Y[j-1])
            score[i][j] = max(scoreSub, scoreDel, scoreIns)
    if toPrint:
        print(X,Y)
        print("Score")
        for row in score:
            print(row)
    lastLine = score[-1]
    return lastLine

# Hirschberg
def Hirschberg(X, Y):
    Z = ""
    W = ""
    # If sequence X has no more char in the seq
    if len(X) == 0:
        for i in range(len(Y)):
            # Concatenate a gap(?) to Z and Y[i] to W
            Z = Z + '-'
            W = W + Y[i]
    # If sequence Y has no more char in the seq
    elif len(Y) == 0:
        for i in range(len(X)):
            Z = Z + X[i]
            W = W + '-'
    elif len(X) == 1 or len(Y) == 1:
        #If one seq has a single use full global alignment. For now there is a biopython implementation of global alignment

        a = pairwise2.align.globalms(X,Y,2,-1,-2,-2)
        Z = a[0].seqA
        W = a[0].seqB

    else:
        xlen = int(len(X))
        xmid = round(len(X)/2)
        ylen = int(len(Y))
        # Ok so this is the meat of the algorithm. The last few conditionals dealt with the base cases.
        # Scoreleft takes NWScore for the first half of sequence X, and Sequence Y

        ScoreL = NWScore(X[:xmid], Y)
        # I imagine this is to start from the right side.
        ScoreR = NWScore(X[xmid:][::-1], Y[::-1])[::-1]

        ymid = numpy.argmax([x+y for (x,y) in zip(ScoreL,ScoreR)])

        # Recursion. This is divide and conquer. Creating partitions that somehow dynamically reassemble even
        # though portions of it are getting reversed?!
        Z = Z+ Hirschberg(X[:xmid],Y[:ymid])[0]+Hirschberg(X[xmid:xlen],Y[ymid:ylen])[0]
        W = W+Hirschberg(X[:xmid],Y[:ymid])[1] +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])[1]

    return Z,W
seq1 = SeqIO.read("seq1.fasta","fasta")
seq2 = SeqIO.read("seq2.fasta","fasta")

Z,W = Hirschberg(seq1.seq,seq2.seq)
print(Z)
print(W)
score = 0
for i in range(len(Z)):
    if Z[i] == '-' or W[i] == '-':
        score+= insertionPenalty
    elif Z[i]  == W[i]:
        score+=match
    else:
        score+=mismatch

print("Score: ",score)

```

# Answer
> 1 votes
**Comparisons** First of all I would simply say pipe it through `Muscle` and then compare the results.

I forget the algorithm advances in the 'new' (now quite old) generation of aligners, but it is worth reading up on the latest Clustal omega, possibly MAFFT and how they view `Muscle`. I do remember is 2018 Clustal Omega upgrade eulogising the historic innovation of Muscle around accuracy and describing how this was implemented. If that works for you, you simply replicate it in your code.

The modern development in aligning algorithms is how they compare in speed versus number of sequences, but accuracy is always measured.

**Amino acids -\> nucleotides** The reason I don't follow the 'accuracy' issues is because its easier to align via amino acids and then switch back to codons. If you are sticking with this gap, penalty system using amino acids is a good approach and then simply back-translating, but if its a UTR then that will not work.

**Brute-forcing** the gap penalties is certainly an option (which is the advantages of coding it). Proving why it has worked might be a bit tricky, but its doable. You could simply fire parameters onto the data and use regex (`import re`) or straight character matching `==` to assess the improvements. Regex is better to assess improvements, whilst `==` approach is much faster. It depends on the parameter space (I think regex). Moving to an coded 'optical analysis' is an easy modification and could be useful.

---

The implication from the comment is there's suboptimal code ...

* The raw default output to Biopython's pairwise should be used to debug this. Then parameters changed stepwise to understand the anomalous output.

Skimming the code I don't spot anything obvious. What I don't understand is why the code has been duplicated

```
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

```

This is present twice. I think Python would complain about the duplication of functions. The duplication of variables is fine and probably imports.

---

A looked through the code a final time and yes there are bugs

```
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

```

That can't be right because whatever you call the functions (which you do), each time the output is `-3`

No matter what you input the output is always `-3`. If this was OOP and these were class methods in some sort of method chain that might be different. These are functions which are operating not different to calling the variable direct. Thus if you substitute

```
... + dele("a")

```

with

```
... + deletePenalty

```

or even

```
.... + in("a") # I know this isn't supposed to be here but its to make a point

```

The answer in all cases is `-3`.

If you asking for bugs, thats got to be involved.

Out of interest are you using an IDE with breakpoints etc...?

# Answer
> 0 votes
Large sequences of dissimilar size? Sounds like you want a local alignment algorithm. Hirschberg is a global alignment algorithm, which will often perform unexpectedly on large sequences of dissimilar size.

---
Tags: sequence-alignment, phylogenetics, phylogeny
---


title: Galaxy vs Cancer Genomics Cloud

text:
# Question
Title: Galaxy vs Cancer Genomics Cloud
I want to use cloud computing services to analyze my oncological RNA-seq data.

I came across **Galaxy** and **Cancer Genomics Cloud** as two cloud computing platforms which are used by the research community.

Can someone give me the pros and cons of using either **Galaxy** or **Cancer Genomics Cloud**.

# Answer
OP update reveals the data is RNA-seq transcript abundance, that is very different from identifying novel driver mutations. So the stuff below isn't relevant.

---

Just to point out the Broad Institute has been using Google Cloud Platform (GCP) for years and they must be amongst the leaders in processing next gen sequencing data on cancer oncogenesis. They do have specific GCP pipelines available via Google.

I've been tracking their work for sometime because of potential 'cross-over'. It is kinda impressive, albeit there are a few things they (excessively) over complicate.

Just to throw it out there, and it depends what you mean by 'data', if I was processing 'driver mutations' in NGS there would be no question I'd use GCP under the Broad Institutes pipeline (they've changed its name, but you can just email them). So just to labor (labour) the point it depends whatever you mean by 'data'
> 0 votes

---
Tags: cloud
---


# Question
Title: How to convert a BAM file into a VCF file
I know that there is a lot of methods to do so, but most of them need to have a reference file. (For example, the GATK force me to provide a reference.) However, I think the BAM file I am working with already embedded the reference information in it. I tried to use the Samtool (mpileup) to do that, but it seems that the vcf file I get has different formats compared with the VCF documentations. Could anyone give me some advice?

Also, could anyone provide a method to validate the vcf file without using a reference file?

Thanks a lot! The following is what I got when running:

```
samtools view -H

```

```
@HD VN:1.5  GO:none SO:coordinate
@SQ SN:chr1 LN:248956422
@SQ SN:chr2 LN:242193529
@SQ SN:chr3 LN:198295559
@SQ SN:chr4 LN:190214555
@SQ SN:chr5 LN:181538259
@SQ SN:chr6 LN:170805979
@SQ SN:chr7 LN:159345973
@SQ SN:chr8 LN:145138636
@SQ SN:chr9 LN:138394717
@SQ SN:chr10    LN:133797422
@SQ SN:chr11    LN:135086622
@SQ SN:chr12    LN:133275309
@SQ SN:chr13    LN:114364328
@SQ SN:chr14    LN:107043718
@SQ SN:chr15    LN:101991189
@SQ SN:chr16    LN:90338345
@SQ SN:chr17    LN:83257441
@SQ SN:chr18    LN:80373285
@SQ SN:chr19    LN:58617616
@SQ SN:chr20    LN:64444167
@SQ SN:chr21    LN:46709983
@SQ SN:chr22    LN:50818468
@SQ SN:chrX LN:156040895
@SQ SN:chrY LN:57227415
@SQ SN:chrM LN:16569
@SQ SN:chr1_KI270706v1_random   LN:175055
@SQ SN:chr1_KI270707v1_random   LN:32032
@SQ SN:chr1_KI270708v1_random   LN:127682
@SQ SN:chr1_KI270709v1_random   LN:66860
@SQ SN:chr1_KI270710v1_random   LN:40176
@SQ SN:chr1_KI270711v1_random   LN:42210
@SQ SN:chr1_KI270712v1_random   LN:176043
@SQ SN:chr1_KI270713v1_random   LN:40745
@SQ SN:chr1_KI270714v1_random   LN:41717
@SQ SN:chr2_KI270715v1_random   LN:161471
@SQ SN:chr2_KI270716v1_random   LN:153799
@SQ SN:chr3_GL000221v1_random   LN:155397
@SQ SN:chr4_GL000008v2_random   LN:209709
@SQ SN:chr5_GL000208v1_random   LN:92689
@SQ SN:chr9_KI270717v1_random   LN:40062
@SQ SN:chr9_KI270718v1_random   LN:38054
@SQ SN:chr9_KI270719v1_random   LN:176845
@SQ SN:chr9_KI270720v1_random   LN:39050
@SQ SN:chr11_KI270721v1_random  LN:100316
@SQ SN:chr14_GL000009v2_random  LN:201709
@SQ SN:chr14_GL000225v1_random  LN:211173
@SQ SN:chr14_KI270722v1_random  LN:194050
@SQ SN:chr14_GL000194v1_random  LN:191469
@SQ SN:chr14_KI270723v1_random  LN:38115
@SQ SN:chr14_KI270724v1_random  LN:39555
@SQ SN:chr14_KI270725v1_random  LN:172810
@SQ SN:chr14_KI270726v1_random  LN:43739
@SQ SN:chr15_KI270727v1_random  LN:448248
@SQ SN:chr16_KI270728v1_random  LN:1872759
@SQ SN:chr17_GL000205v2_random  LN:185591
@SQ SN:chr17_KI270729v1_random  LN:280839
@SQ SN:chr17_KI270730v1_random  LN:112551
@SQ SN:chr22_KI270731v1_random  LN:150754
@SQ SN:chr22_KI270732v1_random  LN:41543
@SQ SN:chr22_KI270733v1_random  LN:179772
@SQ SN:chr22_KI270734v1_random  LN:165050
@SQ SN:chr22_KI270735v1_random  LN:42811
@SQ SN:chr22_KI270736v1_random  LN:181920
@SQ SN:chr22_KI270737v1_random  LN:103838
@SQ SN:chr22_KI270738v1_random  LN:99375
@SQ SN:chr22_KI270739v1_random  LN:73985
@SQ SN:chrY_KI270740v1_random   LN:37240
@SQ SN:chrUn_KI270302v1 LN:2274
@SQ SN:chrUn_KI270304v1 LN:2165
@SQ SN:chrUn_KI270303v1 LN:1942
@SQ SN:chrUn_KI270305v1 LN:1472
@SQ SN:chrUn_KI270322v1 LN:21476
@SQ SN:chrUn_KI270320v1 LN:4416
@SQ SN:chrUn_KI270310v1 LN:1201
@SQ SN:chrUn_KI270316v1 LN:1444
@SQ SN:chrUn_KI270315v1 LN:2276
@SQ SN:chrUn_KI270312v1 LN:998
@SQ SN:chrUn_KI270311v1 LN:12399
@SQ SN:chrUn_KI270317v1 LN:37690
@SQ SN:chrUn_KI270412v1 LN:1179
@SQ SN:chrUn_KI270411v1 LN:2646
@SQ SN:chrUn_KI270414v1 LN:2489
@SQ SN:chrUn_KI270419v1 LN:1029
@SQ SN:chrUn_KI270418v1 LN:2145
@SQ SN:chrUn_KI270420v1 LN:2321
@SQ SN:chrUn_KI270424v1 LN:2140
@SQ SN:chrUn_KI270417v1 LN:2043
@SQ SN:chrUn_KI270422v1 LN:1445
@SQ SN:chrUn_KI270423v1 LN:981
@SQ SN:chrUn_KI270425v1 LN:1884
@SQ SN:chrUn_KI270429v1 LN:1361
@SQ SN:chrUn_KI270442v1 LN:392061
@SQ SN:chrUn_KI270466v1 LN:1233
@SQ SN:chrUn_KI270465v1 LN:1774
@SQ SN:chrUn_KI270467v1 LN:3920
@SQ SN:chrUn_KI270435v1 LN:92983
@SQ SN:chrUn_KI270438v1 LN:112505
@SQ SN:chrUn_KI270468v1 LN:4055
@SQ SN:chrUn_KI270510v1 LN:2415
@SQ SN:chrUn_KI270509v1 LN:2318
@SQ SN:chrUn_KI270518v1 LN:2186
@SQ SN:chrUn_KI270508v1 LN:1951
@SQ SN:chrUn_KI270516v1 LN:1300
@SQ SN:chrUn_KI270512v1 LN:22689
@SQ SN:chrUn_KI270519v1 LN:138126
@SQ SN:chrUn_KI270522v1 LN:5674
@SQ SN:chrUn_KI270511v1 LN:8127
@SQ SN:chrUn_KI270515v1 LN:6361
@SQ SN:chrUn_KI270507v1 LN:5353
@SQ SN:chrUn_KI270517v1 LN:3253
@SQ SN:chrUn_KI270529v1 LN:1899
@SQ SN:chrUn_KI270528v1 LN:2983
@SQ SN:chrUn_KI270530v1 LN:2168
@SQ SN:chrUn_KI270539v1 LN:993
@SQ SN:chrUn_KI270538v1 LN:91309
@SQ SN:chrUn_KI270544v1 LN:1202
@SQ SN:chrUn_KI270548v1 LN:1599
@SQ SN:chrUn_KI270583v1 LN:1400
@SQ SN:chrUn_KI270587v1 LN:2969
@SQ SN:chrUn_KI270580v1 LN:1553
@SQ SN:chrUn_KI270581v1 LN:7046
@SQ SN:chrUn_KI270579v1 LN:31033
@SQ SN:chrUn_KI270589v1 LN:44474
@SQ SN:chrUn_KI270590v1 LN:4685
@SQ SN:chrUn_KI270584v1 LN:4513
@SQ SN:chrUn_KI270582v1 LN:6504
@SQ SN:chrUn_KI270588v1 LN:6158
@SQ SN:chrUn_KI270593v1 LN:3041
@SQ SN:chrUn_KI270591v1 LN:5796
@SQ SN:chrUn_KI270330v1 LN:1652
@SQ SN:chrUn_KI270329v1 LN:1040
@SQ SN:chrUn_KI270334v1 LN:1368
@SQ SN:chrUn_KI270333v1 LN:2699
@SQ SN:chrUn_KI270335v1 LN:1048
@SQ SN:chrUn_KI270338v1 LN:1428
@SQ SN:chrUn_KI270340v1 LN:1428
@SQ SN:chrUn_KI270336v1 LN:1026
@SQ SN:chrUn_KI270337v1 LN:1121
@SQ SN:chrUn_KI270363v1 LN:1803
@SQ SN:chrUn_KI270364v1 LN:2855
@SQ SN:chrUn_KI270362v1 LN:3530
@SQ SN:chrUn_KI270366v1 LN:8320
@SQ SN:chrUn_KI270378v1 LN:1048
@SQ SN:chrUn_KI270379v1 LN:1045
@SQ SN:chrUn_KI270389v1 LN:1298
@SQ SN:chrUn_KI270390v1 LN:2387
@SQ SN:chrUn_KI270387v1 LN:1537
@SQ SN:chrUn_KI270395v1 LN:1143
@SQ SN:chrUn_KI270396v1 LN:1880
@SQ SN:chrUn_KI270388v1 LN:1216
@SQ SN:chrUn_KI270394v1 LN:970
@SQ SN:chrUn_KI270386v1 LN:1788
@SQ SN:chrUn_KI270391v1 LN:1484
@SQ SN:chrUn_KI270383v1 LN:1750
@SQ SN:chrUn_KI270393v1 LN:1308
@SQ SN:chrUn_KI270384v1 LN:1658
@SQ SN:chrUn_KI270392v1 LN:971
@SQ SN:chrUn_KI270381v1 LN:1930
@SQ SN:chrUn_KI270385v1 LN:990
@SQ SN:chrUn_KI270382v1 LN:4215
@SQ SN:chrUn_KI270376v1 LN:1136
@SQ SN:chrUn_KI270374v1 LN:2656
@SQ SN:chrUn_KI270372v1 LN:1650
@SQ SN:chrUn_KI270373v1 LN:1451
@SQ SN:chrUn_KI270375v1 LN:2378
@SQ SN:chrUn_KI270371v1 LN:2805
@SQ SN:chrUn_KI270448v1 LN:7992
@SQ SN:chrUn_KI270521v1 LN:7642
@SQ SN:chrUn_GL000195v1 LN:182896
@SQ SN:chrUn_GL000219v1 LN:179198
@SQ SN:chrUn_GL000220v1 LN:161802
@SQ SN:chrUn_GL000224v1 LN:179693
@SQ SN:chrUn_KI270741v1 LN:157432
@SQ SN:chrUn_GL000226v1 LN:15008
@SQ SN:chrUn_GL000213v1 LN:164239
@SQ SN:chrUn_KI270743v1 LN:210658
@SQ SN:chrUn_KI270744v1 LN:168472
@SQ SN:chrUn_KI270745v1 LN:41891
@SQ SN:chrUn_KI270746v1 LN:66486
@SQ SN:chrUn_KI270747v1 LN:198735
@SQ SN:chrUn_KI270748v1 LN:93321
@SQ SN:chrUn_KI270749v1 LN:158759
@SQ SN:chrUn_KI270750v1 LN:148850
@SQ SN:chrUn_KI270751v1 LN:150742
@SQ SN:chrUn_KI270752v1 LN:27745
@SQ SN:chrUn_KI270753v1 LN:62944
@SQ SN:chrUn_KI270754v1 LN:40191
@SQ SN:chrUn_KI270755v1 LN:36723
@SQ SN:chrUn_KI270756v1 LN:79590
@SQ SN:chrUn_KI270757v1 LN:71251
@SQ SN:chrUn_GL000214v1 LN:137718
@SQ SN:chrUn_KI270742v1 LN:186739
@SQ SN:chrUn_GL000216v2 LN:176608
@SQ SN:chrUn_GL000218v1 LN:161147
@SQ SN:chrEBV   LN:171823
@RG ID:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1  CN:Novogene LB:Wit_Exon9_DEL    PL:illuminPU:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1    SM:Wit_Exon9_DEL
@PG ID:bwa  PN:bwa  CL:bwa mem -t 8 -M -R @RG\tID:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1\tSM:Wit_Exon9_DEL\tLB:Wit_Exon9_DEL\tPU:Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1\tPL:illumina\tCN:Novogene /PUBLIC/source/HW/Disease/database/genome/human/hg38/bwa_index_whole/hg38.fa /TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/QC/Wit_Exon9_DEL/Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1_1.clean.fq.gz /TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/QC/Wit_Exon9_DEL/Wit_Exon9_DEL_FKDN202385493-1A_H5HHJDSXY_L1_2.clean.fq.gz  VN:0.7.17-r1188
@PG ID:sambamba CL:view -S -f bam /dev/stdin    PP:bwa  VN:1.0
@PG ID:MarkDuplicates   VN:2.18.9-SNAPSHOT  CL:MarkDuplicates INPUT=[/TJPROJ6/GB_HUMAN/PJ_AI/service/cromwell/cromwell-executions/ToSomatic/fddd223f-25fd-4533-ae86-14e106c58ec5/call-Mapping/shard-0/Mapping/14f30bea-add3-4ead-b7e1-fcb202ae18f5/call-picardMarkdup/inputs/-714370009/Wit_Exon9_DEL.merge.bam] OUTPUT=/TJPROJ6/GB_HUMAN/PJ_AI/cancer/auto-WES-Somatic/1000/X202SC2005/X202SC20050050-Z01/F002/Mapping/Wit_Exon9_DEL/Wit_Exon9_DEL.final.bam METRICS_FILE=metric.txt REMOVE_DUPLICATES=false ASSUME_SORTED=true OPTICAL_DUPLICATE_PIXEL_DISTANCE=2500 TMP_DIR=[.] VALIDATION_STRINGENCY=SILENT MAX_RECORDS_IN_RAM=4000000 CREATE_INDEX=true    MAX_SEQUENCES_FOR_DISK_READ_ENDS_MAP=50000 MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=8000 SORTING_COLLECTION_SIZE_RATIO=0.25 TAG_DUPLICATE_SET_MEMBERS=false REMOVE_SEQUENCING_DUPLICATES=false TAGGING_POLICY=DontTag CLEAR_DT=true ADD_PG_TAG_TO_READS=true DUPLICATE_SCORING_STRATEGY=SUM_OF_BASE_QUALITIES PROGRAM_RECORD_ID=MarkDuplicates PROGRAM_GROUP_NAME=MarkDuplicates READ_NAME_REGEX=<optimized capture of last three ':' separated fields as numeric values> MAX_OPTICAL_DUPLICATE_SET_SIZE=300000 VERBOSITY=INFO QUIET=false COMPRESSION_LEVEL=5 CREATE_MD5_FILE=false GA4GH_CLIENT_SECRETS=client_secrets.json USE_JDK_DEFLATER=false USE_JDK_INFLATER=false   PN:MarkDuplicates
@PG ID:samtools PN:samtools PP:sambamba VN:1.10 CL:samtools view -H /Users/XXXX/Desktop/data/02.Bam/XXXX.final.bam
@PG ID:samtools.1   PN:samtools PP:MarkDuplicates   VN:1.10 CL:samtools view -H /Users/XXXX/Desktop/data/02.Bam/XXXX.final.bam

```

The output file I get when using Samtools looks like:

# Answer
To go from `.bam` to `.vcf`, it is fairly easy using `bcftools call`:

```
bcftools mpileup -Ob -o <study.bcf> -f <ref.fa> <sample1.bam>
bcftools call -vmO z -o <study.vcf.gz> <study.bcf>

```

Taken from htslib website.
> 4 votes

# Answer
Recently, I also had a similar question and got my answer from Biostars. See my post and the answer I received: Germline variant calling from BAM without reference.

Basically, I was also wondering if it is possible to call SNVs/indels from **aligned** BAM files without reference because reference sequence was already used to produce such files. To this end, it appears some BAM files contain useful (but optional) read-level tags such as `MD` which, if present, would allow one to even reconstruct the entire reference sequence. For more details, please visit the link above.
> -1 votes

---
Tags: sam, vcf, gatk
---


title: Plotting a gene in Seurat

text:
# Question
Title: Plotting a gene in Seurat
I saw in the extensive Seurat documentation for Dimplot (dimensional reduction plot), here, you can plot a gene by specifying it with group.by = "gene" but this does not work in practice.

```
> p <- DimPlot(DG.combined, reduction = "umap", group.by = "Fos")
Error: Cannot find 'Fos' in this Seurat object
> p <- DimPlot(DG.combined, reduction = "umap", group.by = "Inhba")
Error: Cannot find 'Inhba' in this Seurat object

```

Am I mistaken? Do I need to add a metadata column containing expression of the gene? How can I access individual cells gene expression data in seurat so I can create a metadata label based on binned expression levels? Also how would I tell seurat to use a color gradient to label them?

References https://satijalab.org/seurat/reference/dimplot

# Answer
As @swbarnes2 has suggested, you can plot gene expression over the top of a dimensional reduction using the `FeaturePlot` function, using the `features` parameter for gene names:

```
p <- FeaturePlot(DG.combined, reduction = "umap", features = c("Fos", "Inhba"))

```

See the 3k tutorial here:

https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html#cb49
> 2 votes

# Answer
I am not totally sure what you are trying to do, but I think you want FeaturePlot, not DimPlot
> 2 votes

---
Tags: r, rna-seq, scrnaseq, seurat, visualization
---


# Question
Title: tblastn not providing full alignment even on exact-matching sequences
I am using tblastn to search a polypeptide seq against transcripts. One of those transcripts contains the residues for that exact polypeptide once translated, yet BLAST only gets 88% of it rather than it being full length.

Here's the peptide:

```
>Syntaxin
MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDD
KMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQID
YRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRD
LHDMFMDMAMLVESQGEMIDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGV
IGGTLG

```

And the transcript.

Here's the first-frame translation of the transcript which contains the exact polypeptide above:

frame +1

```
> - 669
PATRNITIVSIYDYH*NNDLKCVFTE*MEILNRIWILLNSHSP*AQRLHAI**MTKDRLAALKAAQSDDDD
NDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKHSAILSAPQTDDKMKEELEELMSEIKKNAN
KVRAKLKVIEQNIEQEEHTNKSSADLRIRKTQHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGK
TTTNEELEDMIESGNPAIFTQGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEM
IDRIEYNVEQAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG***PSSLPLVYS
LIQLFYAACTTAASNIWTHHNRNSYYYNSNSNSNSNCYSNYFSNYYSNRKHYYYYHYYYYYYYYYYYYYYY
YYYYYYDYYYYYYYYYYYYYYYYYCYYYYYYYYYHFHYYYHHYYYHYYYYHYYYHYHHHHHYHYHYHYHHH
HHYHYHYHYHYHYYYYYYYYYYYYYYYYYYYYYYCYYYYYYTFSLVISLYGSSEM*MRLLVLSLYTFYGLF
YPLQQKGSFFF*SLQQTGTGFYNFFSYCQLHNLSIPSPSYFF*EWVFSFLKKHLIEFFCFLILLKFVLSFS
PIISLIISVCRKKISRIFYTIYLKMS**NHYLLF

```

Yet when I do a local tblastn (using the modern blast+) or use NCBI's pairwise blast tool it doesn't return a full-length match.

Why does it leave off the last 12% of the sequence match?

Here's the alignment:

# Answer
That's because by default, blast masks low complexity regions. With the online blast2seqs tool, you can disable this by clicking on the "Algorithm options" button:

And then deselect the "Low complexity regions" filter:

This will let you align the entire sequence as expected:

```
Alignment statistics for match #1
Score   Expect  Method  Identities  Positives   Gaps
Frame
588 bits(1516)  0.0 Compositional matrix adjust.    290/290(100%)   290/290(100%)   0/290(0%)
+1
Query  1     MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  60
             MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH
Sbjct  154   MTKDRLAALKAAQSDDDDNDDVAVTVDSSGFMEEFFEQVDEIREMIDKIASNVDEVKKKH  333

Query  61    SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  120
             SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT
Sbjct  334   SAILSAPQTDDKMKEELEELMSEIKKNANKVRAKLKVIEQNIEQEEHTNKSSADLRIRKT  513

Query  121   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  180
             QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT
Sbjct  514   QHATLSRKFVEVMNDYNACQIDYRERCKGRIKRQLAITGKTTTNEELEDMIESGNPAIFT  693

Query  181   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  240
             QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE
Sbjct  694   QGIIMETQQANETLADIEARHNDIMKLETSIRDLHDMFMDMAMLVESQGEMIDRIEYNVE  873

Query  241   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  290
             QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG
Sbjct  874   QAVDYIETAKMDTKKAVKYQSKARRKKIMILVCLAILIIILVGVIGGTLG  1023

```

This is the part of your transcript sequence that was masked:

```
acgctaggctgaccctcctccctccctctagtctactccttaatacaattgttctacgccgcgtgcacc

```

I am guessing it is the repeated `ccct` nucleotides that resulted in the sequence being masked, but I am not sure.
> 2 votes

---
Tags: blast
---


# Question
Title: Nextflow: dynamically setting resource directives for a process
I'm trying to test variations of CPU and memory for a specific process (`samtools sort`) to measure computational efficiency (measured only by time). To do this, I'd like to set the CPU and memory for a `Nextflow` `process` dynamically.

**Ultimate question:** How can I dynamically set CPU and memory for many instances of a single `Nextflow` process definition?

# Attempt 1

The most obvious solution was to submit distinct `Nextflow` workflows for each combination of CPU and memory by specifying the CPU and memory on the command line (i.e., when calling `nextflow run ...`), where my `Nextflow` workflow has global parameters `params.samtools_total_threads` and `params.samtools_mem_per_thread`.

This didn't work perfectly because there appeared to be 'collisions' between parent `Nextflow` instances trying to lock the history file. I got the following error:

```
Can't lock file: /path/to/project/.nextflow/history -- Nextflow needs to run in a file system that supports file locks

```

I could try delaying submissions, but I suspect that will still result in collisions.

Here's the `SLURM` script that I'm using:

```
TOTAL_THREADS=$1
MEM_PER_THREAD=$2

nextflow run 02-SAMTOOLS_SORT.nf \
    -with-report "reports/report.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.html" \
    -with-trace "traces/trace.${TOTAL_THREADS}_threads.${MEM_PER_THREAD}_GB_per_thread.txt" \
    --samtools_total_threads $TOTAL_THREADS \
    --samtools_mem_per_thread $MEM_PER_THREAD

```

My workflow had the following global parameters:

```
/*
 * Make this pipeline a nextflow 2 implementation
 */
nextflow.enable.dsl=2

...

params.samtools_total_threads = '4'
params.samtools_mem_per_thread = '8'

...

```

# Attempt 2

After Attempt 1 failed, I figured the most 'Nextflow-like way' would be to generate tuples for each combination of CPU and memory that I want to test. Unfortunately, it doesn't appear that it's possible to set CPU and memory directives using parameters passed into the `process`, but only with global variables.

I get the following `No such variable` error with the below code (note: most of this code is courtesy of @Steve):

```
N E X T F L O W  ~  version 21.04.3
Launching `02-SAMTOOLS_SORT.nf` [desperate_aryabhata] - revision: fb9e79055a
No such variable: total_threads

 -- Check script '02-SAMTOOLS_SORT.nf' at line: 33 or see '.nextflow.log' file for more details

```

**Process code**

```
process samtools_coordinate_sort_proc {

    tag { "${bam.baseName}" }

    executor='slurm'
    queue = 'normal'
    cpus = total_threads
    memory = total_mem.GB
    clusterOptions = "--time 5:00:00 --account <account>"

    /* delete files upon completion (I think) */
    cleanup = true

    input:
    path(bam)
    val(total_threads)
    val(total_mem)

    output:
    path("${bam.baseName}.csorted.bam")

    script:

    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    println "avail_mem: $avail_mem"
    def mem_per_thread = avail_mem ? "${avail_mem}G" : ''
    println "mem_per_thread: $mem_per_thread"

    """
    samtools sort \\
        -@ "${additional_threads}" \\
    -m ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
    -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

```

# Possible attempt 3

Of course, it would be *super* easy to do this with a couple of `bash` scripts, but I'm trying to do it in `Nextflow` for two reasons:

1. To learn `Nextflow` better
2. Because the traces provided by `Nextflow` are super convenient (though I could get this information from `SLURM`, too).

# Answer
It's possible to specify dynamic computing resources using a closure. For example, this sets the cpus and memory directives using the values supplied via the input declaration:

```
nextflow.enable.dsl=2

process test_samtools {

    echo true

    cpus { cpus }
    memory { memory }

    input:
    tuple path(bam), val(cpus), val(memory)

    script:
    def additional_threads = task.cpus - 1

    def avail_mem = task.memory ? task.memory.toGiga().intdiv(task.cpus) : 0
    def mem_per_thread = avail_mem ? "-m ${avail_mem}G" : ''

    """
    echo samtools sort \\
        -@ "${additional_threads}" \\
        ${mem_per_thread} \\
        -o "${bam.baseName}.csorted.bam" \\
        -T "${bam.baseName}.csorted" \\
        --write-index \\
        "${bam}"
    """
}

workflow {

    cpu_ch = Channel.of( 1, 2, 3, 4, 5)
    memory_ch = Channel.of( 1.GB, 2.GB, 4.GB )

    Channel.fromPath( './data/*.bam' )
        | combine( cpu_ch )
        | combine( memory_ch )
        | test_samtools
}

```
> 6 votes

---
Tags: nextflow
---


title: use Kallisto in galaxy

text:
# Question
Title: use Kallisto in galaxy
I want to use Kallisto for sequence alignment in Galaxy. Its description is:

> a program for quantifying abundances of transcripts from bulk and single-cell RNA-Seq data, or more generally of target sequences using high-throughput sequencing reads.

However, I found this field empty:

> **there is no option available for the reference transcriptome.**

* How I use a reference sequence?
* Can I leave the other configuration by default?

# Answer
You have to click on the dropdown menu and select "Use transcriptome from history", which will allow you to specify your own file.

I suspect the default will be built-in human transcriptome.

Solution is not tested, so let us know if it actually works.
> 1 votes

---
Tags: sequence-alignment, transcriptome, galaxy, kallisto
---


# Question
Title: Programmatic way of accessing the FASTA sequence similarity tool (or similar) in Python
I am looking for a tool that performs a sequence similarity algorithm for proteins. More specifically, I am looking for something that would be usable in Python (or anything else usable in the command line actually) which has a similar output to the FASTA tool, i.e:

For my use case, the 'Source' column is not super relevant, so it does not matter if it's not included in the output. What's most important is the DB:ID and the match statistics.

Have you used anything similar? Do you have any suggestions regarding such a tool? Thanks in advance!

# Answer
> 2 votes
Apparently there is a REST API for FASTA. Looks like there are few dependencies but other than that should be straightforward to use.

---
Tags: sequence-alignment, fasta, python, proteins
---


# Question
Title: Querying The Cancer Genome Atlas (TCGA) for gene expression
The UCI gene expression cancer RNA-Seq dataset (https://archive.ics.uci.edu/ml/datasets/gene+expression+cancer+RNA-Seq) consists of 5 types of cancers (BRCA, KIRC, LUAD, COAD, PRAD) and gene expression levels of 20,531 genes for 801 samples. The dataset looks like this:

How to query in The Cancer Genome Atlas (TCGA) to create an updated version of this dataset?

I've done the following query.

Is the given query and the UCI dataset equivalent? If not, how to query properly?

P.S.: I'm a CS student; do not have any biology background.

# Answer
To start with, your query above is selecting for breast cancer samples (TCGA-BRCA) only. The original dataset seems to contain 4 other cancer types.

The link you provided has a link to https://www.synapse.org/#!Synapse:syn4301332 and by looking at filenames of the 14 files there, I can tell that there are different types of data: miRNA expression, gene expression, methylation, ... And your query is limited to gene expression only. I am not even sure if you can use just a single query to fetch all of the 14 data files, you might need to use multiple queries.

If you are familiar with R, I would recommend TCGABiolinks, it will abstract away a lot of otherwise annoying manual steps. Apparently you can also use Python retrieve TCGA data:https://docs.gdc.cancer.gov/API/Users\_Guide/Python\_Examples/
> 2 votes

---
Tags: genome, database, gene-expression, cancer
---


title: Understanding fast5 format

text:
# Question
Title: Understanding fast5 format
I would like to reproduce the results on figure 1 from THIS publication. I downloaded the data and explore it with h5ls and h5dump based on HDF5. I would like to understand:

1. `h5ls -r <_fast5file_> | grep Dataset`:

```
/Analyses/Basecall_2D_000/BaseCalled_2D/Alignment Dataset {6925}
/Analyses/Basecall_2D_000/BaseCalled_2D/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Events Dataset {5576}
/Analyses/Basecall_2D_000/BaseCalled_complement/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_complement/Model Dataset {4096}
/Analyses/Basecall_2D_000/BaseCalled_template/Events Dataset {5793}
/Analyses/Basecall_2D_000/BaseCalled_template/Fastq Dataset {SCALAR}
/Analyses/Basecall_2D_000/BaseCalled_template/Model Dataset {4096}
/Analyses/Basecall_2D_000/HairpinAlign/Alignment Dataset {5107}
/Analyses/Basecall_2D_000/Log Dataset {SCALAR}
/Analyses/Calibration_Strand_000/Log Dataset {SCALAR}
/Analyses/EventDetection_000/Reads/Read_0/Events Dataset {11465}

```

Alignment? Then `complement/Fastq` is not the complement of `template/Fastq` (analyzed with h5dump). What is it? Also, events have sometimes 4

2. Some other datasets that were acquired somewhere else have the raw data and fastq but there is no associated information about kmer and mean(current)/std/length/etc... How do I generate events from raw data?

I went through all documentations and tutorials that I could found in the web but couldn't find anything simple and complete. It feels as if it was a secret.

Example *fast5file* would be `ecoli_er2925.MSssI.timp.100215.fast5/pass/kelvin_100215_ecolisssI_2515_1_ch195_file0_strand.fast5` from PRJEB13021

# Answer
FAST5 is a proprietary format developed by Oxford Nanopore Technologies, and there's not much great documentation online. All FAST5 files will have the `Raw/` field, which contains the original measured `Signal`. Additional `Analyses/` fields can be added by tools such as basecallers which convert signal-\>ACGT (Guppy) or signal-level aligners (Tombo).

In particular, `/Analyses/Basecall_2D_000/` contains the results of 2D basecalling, added by Guppy when run with the `--post-out` flag. I think Tombo adds the event-level annotations you mention (k-mer mean/std/length). In 2D basecalling, both the physical forward and reverse strands for the same stretch of DNA are measured and basecalled. The basecalls will contain errors, so these will not be exact reverse complements of one another.

Hope this answers your questions, I have a slightly more in-depth explanation on my blog: https://timd.one/blog/genomics/file\_formats.php
> 2 votes

# Answer
I'll give a different perspective. Oxford Nanopore FAST5 is NGS within a `HDF5` database. A clear description is here.

HDF5 is a format used for rapid I/O of large data sets.

Oxford Nano appear to have written a bit boilerplate around HDF5 and called it something that sounds like 'fasta'.

A cool part of HDF5 (not FAST5) is `HDF5 View`, which can be downloaded here. Its cool because you can visualise the file structure of HDF5 db via a GUI. I very much suspect if you use HDF5 View on your FAST5 database you will be able to navigate (see not change) the entire file structure.

To answer your question directly `h5ls` is listing the HDF5 database. To construct a HDF5 db you assign the data under a series of 'keys' which forms the directory structure. However, you perhaps might not need Oxford tools, because Python pandas have direct series of methods to connect, construct and manipulate HDF5.

I raised the HDF5 in a question then answered here on Bioinfo SE here.

Again the reason for HDF5's importance is unlike an SQL connection, data transfer and disconnect, HDF5 shifts data very quickly.

I do question how Oxford Nanopore can lay claim to FAST5 because HDF5 is open source and always has been. If they want to lay claim to some boilerplate, be my guest - I write the stuff all the time ;-)

HDF5 however is a serious project with major development time and appears to still be under active development.
> 1 votes

---
Tags: nanopore, fast5, hdf5
---


# Question
Title: Samtools sort: most efficient memory and thread settings for many samples on a cluster
We're preparing to analyze thousands of `.bam` files beginning with re-alignment, sorting, etc.

**Complicated question:** Has anyone investigated the optimal thread and memory settings to sort hundreds to thousands of `.bam` files using `samtools sort` on a cluster?

By 'optimal', I'm looking to sort as many samples as possible in the shortest time frame. I'm mostly interested in optimizing `samtools sort` because the other steps (e.g., alignment) are much more obvious to fine-tune.

For a single sample, throwing large amounts of memory + threads is trivial, but that isn't scalable for thousands of samples. There is also likely a diminishing return at a certain point, even for a single sample.

Ideally, I'd like to see plots identifying the 'sweet spot' for threads + memory for a single sample (e.g., x == threads, y == memory, z == time), and then plots maximizing throughput for more samples. I would expect the data to follow a -log function, where we could identify the optimal settings at the knee in the curve.

There are obviously many other factors (e.g., cluster size, administrative scheduler settings, etc.), but I suspect findings on a given cluster will generalize to other clusters.

Really appreciate any feedback.

# Answer
> 6 votes
# Update 2

Ok, for what I expect will be my final update for this answer, I compared the following:

1. `samtools 1.14` \+ `zlib`
2. `samtools 1.15` \+ `libdeflate` (this is a different version, but shouldn't have a major effect. They released a new version just before I asked cluster admins to compile against `libdeflate` and didn't notice until later.)
3. `sambamba v0.8.2`

**Note**: This time I used the following settings:

1. Only two samples instead of three
2. Max `Nextflow` queue size of 30 to avoid too many threads reading from the same two files.
3. I tested CPU requests from 1 to 7 (step size 1) and then from 9 to 17 (step size 2).
4. For jobs where I allotted few CPUs and low memory, I provided a 200% buffer to prevent the jobs from failing with `OUT OF MEMORY` errors. i.e., I gave the `SLURM` job 2x the amount of memory than I told the tools they could use. `sambamba` was the worst offender for exceeding memory in these conditions, but I gave them all 200% buffer to keep comparisons equal.

Seems pretty clear to me that `sambamba` is the fastest and uses requested CPU resources most efficiently (i.e., is able to use all of the CPUs requested up until ~8 or 9 CPUs). For reference, **the difference between `sambamba` and `samtools + libdeflate` with 9 CPUs is ~10 minutes on average.**

We will certainly be using `sambamba` going forward, probably with ~9 CPUs and ~9GB per thread since I had to provide a 200% buffer in lower CPU/Mem jobs to avoid job failure, anyway.

Again, all code and results are available on my GitHub repository (Samtools sort optimization test).

## (Angle 1) Comparison between `samtools + zlib`, `samtools + libdeflate` (Lsamtools), and `sambamba`

Here, I'm showing the general trend that `real time` decreases with increased CPUs for all three tools, but Mem Per CPU doesn't have a huge effect. `sambamba` out performs `samtools + zlib` regardless of settings. `samtools + libdeflate` out performs `samtools + zlib` until ~11 CPUs, where they perform the same (plot 3 below explains why).

**Interesting finding:** `samtools + libdeflate` performs better than `sambamba` with a single thread, but `sambamba` takes over from there as CPUs increase (again, plot #3 explains why).

## (Angle 2) Comparison between `samtools + zlib`, `samtools + libdeflate` (Lsamtools), and `sambamba`

Here, I'm showing that `samtools + zlib` and `sambamba` performance flattens out at ~11 CPUs (again, plot 3 below explans why).

## Requested CPUs vs. CPU utilization for `samtools + zlib`, `samtools + libdeflate` (Lsamtools), and `sambamba`

Here, I'm showing how well each tool is able to utilize the CPUs allotted to it. It was fascinating to see that `samtools + libdeflate` flattens out very early in it's CPU utilization (i.e., ~3 CPUS) and struggles to cross the 400% CPU utilization with even \>10 CPUs. This is why `samtools + libdeflate` performance flattens out quickly.

`samtools + zlib` does a much better job utilizing allotted CPUs, but still lags behind `sambamba` throughout (and is slower, anyway).

# Update 1

My first plots (below) were showing `duration` on the `z-axis`, which is the time from job *submission* to *completion* rather than job *start* to *completion*. That's why the larger memory requests were much longer in the original plots. `Duration` is still a useful measure, but I'm now including `real time`, which is from *start* to *completion*.

Judging solely by `real time` (above), there is a clear gain between 5 & 7 CPUs (~10 minutes, on average when compared to the same mem-per-thread).

More data coming, including `samtools` compiled against `libdeflate` and `sambamba`.

---

# Original post

I made time to semi-formally look at this. I'm glad I did because the results surprised me. I don't answer all of my questions directly, but these results are good enough to get me going in the right direction.

My code and results are in my `GitHub` repository (Samtools\_sort\_optimization\_test) where you will also find an interactive plot you can download to see the results more clearly (see `SAMTOOLS_SORT/samtools_sort_CPU_and_memory_comparison.tar.gz` in the repository). It's hard to truly appreciate the 3D `plotly` plot until you start rotating it.

**NOTES:**

1. Even though this is from only three samples, I think my results are accurate. I welcome any feedback in case there's a bug.
2. In my experience, `samtools` often exceeds its allotted memory (at least on a `SLURM` scheduler). i.e., the job gets killed by `SLURM`. To avoid this, I requested an additional 20% of memory to provide a buffer for `samtools`. I don't know if this discrepancy is because of `SLURM` or `samtools`, though.
3. These results are using `samtools` compiled against `zlib`. We are preparing a version of `samtools` compiled against `libdeflate`, per @user172818's suggestion. I will add those results once completed. I may also compare against `sambamba`, per @acvill's suggestion.
4. These results are related to my other question about the best sorting/merging method. I will post those results in that question and provide a brief summary here.

# tl;dr (too long, didn't read)

## Primary conclusions based on the plots below

1. **Optimal CPUs** (results as expected): There's a significant gain with some additional CPUs, but it levels off around 7 CPUs (total, not 'additional'; i.e., `samtools sort -@ 6`). I thought we might still see *some* nice gains above 7, but this is probably what most would expect (e.g., similar to what @user172818 suggested).
2. **Optimal memory per thread** (~~results NOT as expected~~): ~~There is very little to gain with more than 1GB of memory (per thread); that surprised me. There are *some* performance gains up until ~7GB per thread (~2-4 minutes faster). What really surprised me, however, is that performance dramatically **decreases** as memory per thread goes \>7GB~~. The updated plot (above) shows there is some gain with \>1GB mem per thread. When running lots of samples at once (hundreds to 1000), mem per thread seems to play a much larger role (data coming).
3. **Effect from sorting status** (results as expected): Largely in an attempt to do this as quickly as possible (laziness?), I used three samples, where one was from an aligned and unsorted `.bam` file, and the other two were from `.bam` files that had been aligned and sorted, but I then collated them (`samtools collate`), rather than re-aligning them. Unsurprisingly, the collated samples (not the average use case) sorted slightly faster than the unsorted sample (average use case).

## Initial opinion

1. **"Small" sample set**: For a small sample set, I would probably use 7 total CPUs and 7GB per thread.
2. **"Large" sample set**: I will probably use 7 CPUs (maybe 5) with only 1GB memory per thread. That will allow the scheduler to keep as many jobs going as possible.

## `samtools sort`: duration vs. CPUs and mem per thread

As mentioned above, an interactive version of this plot is available on my `GitHub`.

**Angle 1**

**Angle 2**

# Additional details

Analyses were run on University of Kentucky's MCC Cluster, which uses `AMD 7702P` processors, using `samtools v.1.14`

# Answer
> 8 votes
I usually use `samtools sort -@4 -m4g`. Samtools will be faster with more threads and more RAM, but I *feel* the gain is not significant. I haven't done a proper evaluation, but I have run samtools sort with different settings from time to time on clusters when I have to consider speed-resource tradeoffs.

When you use samtools, make sure you compile it with libdeflate, which is about twice as fast as zlib on compression. I don't know if the bioconda version built with libdeflate. You may consider the precompiled binary here.

Samtools sort, especially with libdeflate, is fairly fast these days. Performance metrics in the other answer are outdated.

---
Tags: sequence-alignment, sequence-analysis, samtools, benchmarking
---


# Question
Title: Compare mapped reads from BWA -MEM and STAR from .bam files
I want to find and compare the results from STAR and BWA- MEM mapped. I have 150bp paired end reads in sorted.bam files in each case and i want to find in each case uniquely mapped reads, number of reads mapped to multiple loci and number of unmappped reads.

Especially for uniquely mapped reads i tried different solutions from biostars and this Obtaining uniquely mapped reads from BWA mem alignment but i dont know what to trust,

Could you please suggest the best way to find mapping results from bam paired end files ?

# Answer
> 2 votes
Bam files contain a MAPQ field that informs on whether an alignment is unique or not. With STAR, a MAPQ value of 255 means the read is mapped uniquely and a score lower than 255 means not unique.

You can use samtools on a bam file to count or output the uniquely or multimapped reads by filtering for the MAPQ as described here: https://groups.google.com/forum/#!topic/rna-star/ipHc9FrgaPU

For example, count uniquely mapped reads: `samtools view -c -q 255 myBamFile.bam`
Output uniquely mapped reads: `samtools view -q 255 myBamFile.bam`

There are differences between STAR and BWA MAPQ output though. see here for details: https://sequencing.qcfail.com/articles/mapq-values-are-really-useful-but-their-implementation-is-a-mess/

---
Tags: linux, bwa, star
---


# Question
Title: Questions about microbial 16s rRNA phylogenetic, ANI Taxonomy, and GGDC
I have a few questions about these three methods.

1. We commonly use 16s rRNA to identify species and construct a phylogenetic tree. If there's a new isolate from a new and undiscovered ecosystem, should we confirm this new isolate identity with ANI taxonomy/ genome to genome distance calculation (GGDC)?

* If that isolate has a high similarity with a certain species in the 16s rRNA analysis result, should we proceed to analyze further with ANI/GGDC? Or just let it be because we have a robust 16s rRNA analysis? How can we be sure with our 16s rRNA BLAST result?

2. When should we use 16s rRNA rather than ANI/GGDC and vice versa?
3. Is it always necessary to confirm 16s rRNA analysis with ANI/GGDC? (Assuming there's a lot of closely-related genome data available in databases)
4. Should we confirm the ANI result with GGDC and vice versa?

# Answer
I do understand that you are working with is an isolate and thus what occurs in metagenomics you feel is over doing it, particularly if 16S is e.g. \>99% similarity via Blast.

Before reading further you might check my answer history to understand the angle I am approaching this question from. I would see select given methods as the 'gold-standard' other methods are compared against, in addition to single versus multi-gene trees. There is good theoretical criteria for believing this. For example, bacteria are often subject to heavy %GC bias and this makes other methods error prone, because it introduces compositional bias (below).

**Summary** Blast alone - thats a 'no' (theory section). The problem with any one approach is that errors can arise from many different sources stemming from both biology and methodology. Thus a combination of approaches is required. Therefore your maximum approach is a minimum standard.

I personally would go one step further than your question: full maximum likelihood and Bayesian phylogenetics with bootstraps in addition to Blast (minor), ANI and GGDC.

**History** From tradition 16S is the major standard, but even then multi-gene phylogenies were required. Thus genome to genome distance calculation (GGDC) with ANI is essential at a minimum. 16S is back in fashion because its easy to perform a metagenomics study of species population dynamics.

Species delineation has always been problematic in bacteria and the rule of thumb was 5% divergence defined a species boundary.

**Advice** I would use a multi-locus appproach under a full phylogenetics analysis, where a species boundary is defined by \>5% divergence and bootstraps \>80% with Bayes P values \> 0.95.

Recombination is essentially ignored in taxonomy providing there isn't over reliance on a single gene tree, but using BLAST, ANI, GGDC as a method data-mining would justify a select sample size. The classification is then nailed with a phylogenetics analysis.

One of your methods (probably ANI) might be based machine learning/deep learning (its not GGDC) and the accuracy of misassigning a species under this method is high. Historically it was really high (80% accuracy), this might have improved to 90%. Its not like image recognition where the error rate is beyond human perception.

---

**Theory** A fundamental issue is that recombination in bacteria exists on a spectrum from clonality to panmixia. This makes single gene taxonomy risky, thus at theoretical level would not believe 16S was a singularly criterion, but it never was. Your bacteria is unknown (to us) so we'd have to assume some level of recombination.

Compositional bias is well known for taxonomic misassignment. It is readily demonstrated in simulation. Thus if two unrelated bacteria demonstrate e.g. 75% GC, they are more likely to artifactually grouped together on e.g. a distance criterion.
> 4 votes

---
Tags: phylogenetics, phylogeny, bacteria, 16rrna
---


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

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

# Answer
No matter how you approach this problem 110+ GB is big:

Your options are:

1. Break up the file and perform the intersection on a chunk by chunk basis
2. Ask your system administrator to install the program
3. 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.
> 1 votes

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

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

---
Tags: bedtools
---


title: How to install DnaSP on a Mac

text:
# Question
Title: How to install DnaSP on a Mac
I need to install DnaSp, a nice GUI originally written for Windows for the population genetics of nucleotide sequence data.

However, its not working on my mac. I'm using wine to install it using these instructions for wine. Wine usually successfully ports Windows programs.

it successfully installs but then I get this error

```
Mollies-MacBook-Air:dnasp51001 molliepassacantando$ WINEPREFIX=~/.wine64 wine regsvr32.exe scrrun.dll mfc40.dll threed32.ocx
regsvr32: Successfully registered DLL 'scrrun.dll'
regsvr32: Successfully registered DLL 'mfc40.dll'
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
regsvr32: Successfully registered DLL 'threed32.ocx'
Mollies-MacBook-Air:dnasp51001 molliepassacantando$ winedevice.exe(30828,0x405ed000) malloc: *** error for object 0x401cea08: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug

```

and cannot run it.

Has anyone successfully installed DnaSp on a Mac or know of an alternate program I could use? I need to be able to calculate haplotype frequencies.

# Answer
From the comments @MollePasscantando successfully ported DNAsp onto a Mac (Air) via Wine, but was not sure how. This is good news because Wine is very commonly used for porting Windows exes. The website states,

> Oracle VirtualBox (free open-source). Parallels, CrossOver, VMware Fusion (paid-versions; trial period)

From the website it was update in February 26, 2019 and its current Beta Version is 6.12.03.

There was a major published update to the package in Mol Biol Evol in 2017 (version 6) for large data sets (NGS).

Thus the question is still relevant because they still have not written an OSX specific version and this will limit its performance for the data sets it is now targeted for.
> 0 votes

---
Tags: genome, phylogenetics, software-installation
---


# Question
Title: Are there any "canonical" sources of C. elegans connectome?
I found that there are a couple of projects that allow the user to simulate the neural system of *C elegans*, for example OpenWorm one update being here.

I don't know how to navigate sources in bioinformatics and thus struggle to evaluate their quality.

**Question** Are there any fundamental / most widely accepted data sources for raw data of C Elegans connectome or electron microscope scans of this organism, similar to Human Genome Project? If not, what are the ways to derive the data?

# Answer
> 1 votes
*Answer from @terdon, converted from comment (please improve if possible):*

Have you checked on WormBase? That's the main repository for c. elegans data. I don't know that they have connectome data, but if they exist, that's where they will be stored.

*\[see here for more information on the purposes of comments in the StackExchange network\]*

---
Tags: public-databases
---


# Question
Title: Index VCF file for rapid access by rsID?
With tabix one can index a VCF file for rapid variant retrieval based on genomic position. I'm wondering if there are any tools that will index a VCF file to allow rapid retrieval using rsIDs and/or other metadata? I'm aware of awk/grep/vcftools one-liners for this purpose, but I'd like to avoid scanning a huge VCF each time I need to retrieve the coordinates of a new batch of rsIDs.

# Answer
John Marshall said:

> *Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.*

I just published rsidx, a new package for indexing VCF files by rsID. As John suggests, I simply created a mapping from rsID to genomic coordinates. These are stored in an sqlite3 database.

It currently takes...several...hours to index build 151 of dbSNP on GRCh38. I'm looking into optimizing sqlite3 parameters to accelerate the indexing, as well as an alternative approach based on a two step process (populate an in-memory sqlite3 db, then dump DB to disk). But in the meantime, the search capability works great.
> 0 votes

# Answer
I am not aware of any widespread index formats for query-by-QNAME for BAM or query-by-ID for VCF. I thought there were historical *samtools{-devel,-help}* threads to the effect of “in principle one could index a named-sorted BAM file, but there's insufficient demand so no-one's ever implemented it”, but I can't find any just now.

Moreover note that VCF files are by definition sorted by genomic position, so such an index would just be a hash table of rsID to file offset or genomic position.

If you are using dbSNP, probably you have a dbSNP database or API that you can use to query by rsID and get genomic position back. So you can retrieve your VCF contents by rsID in two stages: do that dbSNP query, and use its results to find the variants in your VCF file by position.
> 2 votes

# Answer
**high-perf-bio.**

I present my *high-perf-bio* project, which, for the most part, solves the problem of fast extraction of anything values from VCF, BED and unformatted data.

**Preparation.**

1/ Run the automatic downloading of the toolkit and installation of the dependencies.

```
wget https://github.com/PlatonB/high-perf-bio/archive/refs/heads/master.zip && unzip -o master.zip && rm master.zip && cd high-perf-bio-master && bash __install_3rd_party.sh

```

2/ Reboot.

3/ `cd path/to/high-perf-bio-master`

**Create DB with rsIDs.**

```
python create.py -S path/to/gzipped_tables_dir -i rsIDs_field_path

```

**Search rsIDs set in DB.**

```
python annotate.py -S path/to/annotated_gzipped_tables_dir -D db_name -T path/to/trg_dir #-c rsIDs_column_num -f rsIDs_field_path

```

P.S.1. Detailed help for each *high-perf-bio* component: `python tool_name.py -h`

P.S.2. For VCF and BED you can specify the minimum of arguments.

P.S.3. Readme in English is in the early stages of writing. If you have any questions, write to Issues.
> 2 votes

# Answer
Take a look at `grabix`. Using `grabix` you could implement a binary search on a VCF sorted by rsID and compressed using `bgzip`.
> 1 votes

---
Tags: vcf, file-formats, variation, indexing
---


title: Definition of Copy Number Abundance

text:
# Question
Title: Definition of Copy Number Abundance
Reading from this paper:

> The remaining 9 new copy-number alterations showed altered coverage patterns on whole-genome sequencing but could not be confirmed by orthogonal methods because of their small size, **low abundance**, or both (Fig. 2C, Fig. S2C, and Table S5)

In the supplementary file, CNA abundance as percentage is defined as:

$$Abundance = \frac{2^{L2R} - 1}{\frac{CN}{2} - 1}\*100$$

where L2R is the log$\_{2}$ normalized coverage ratio vs. a panel of normals and CN is the estimated copy number for the event.

I haven't found a better definition/explanation behind "CNA abundance", is it defined elsewhere, or is it a common term in CNA analysis?

# Answer
> 0 votes
My guess by analyzing the definition is that CNA abundance is a proxy for the percentage of (tumor) cells carrying the alteration.

---
Tags: copy-number
---


# Question
Title: What is the correct order of flooring-normalization-batch correction for microarrays?
*This question was also asked on Biostars*

I am trying to learn and understand the correct order of data processing steps for microarrays.

I have data which already was analyzed by a researcher using Partek (with an ANOVA model with batch effects). I want to test the results versus a specific gene set in GSEA. However, as I was informed, the genes were not filtered for low expression before ANOVA. So I want to filter genes out with low expression as is recommended for GSEA. However, I wonder whether it is legitimate to filter after DE was already done.

In what order should one apply flooring (setting all values below a threshold to that threshold's value), normalization, and batch correction?

# Answer
The most important thing to do first is batch correction, and it looks like that was done in this instance.

The ordering in this case seems reasonable. Normalisation is usually applied to a full dataset (it depends on the assumed model that was used for normalisation, and I'm not familiar with how Partek does things).

Filtering after DE should be fine. I'm at least aware that low-count filtering with DESeq2 \[cDNASeq rather than microarray\] shouldn't influence differential expression results, but it can speed up the calculations slightly:

> While it is not necessary to pre-filter low count genes before running the DESeq2 functions, there are two reasons which make pre-filtering useful: by removing rows in which there are very few reads, we reduce the memory size of the dds data object, and we increase the speed of the transformation and testing functions within DESeq2.

https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#pre-filtering
> 1 votes

---
Tags: normalization, microarray, batch-effects
---


# Question
Title: BWA mem | samtools view: Intermittent parsing error
# Update

The issue was that `bwa` was running out of memory and failing, but that error wasn't floating to the top (see @Steve's answer, below). I was getting an error from `samtools`. I should have been able to figure that out, but hopefully this post saves someone else the trouble.

# Original post

I'm experiencing a strange intermittent parsing error when piping `bwa mem` to `samtools view`. It only seems to happen when the command runs within a `Nextflow` pipeline (though it seems unlikely that `Nextflow` has anything to do with it). If I `cd` to the directory and run `bash .command.sh`, it runs to completion without issue.

# Nextflow process

Here's the `Nextflow` process I'm running (code credit to @Steve, though modified to handle `.bam` or `.cram`). The portion that's running in this case

```
process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    if [ "$params.output_format" = "bam" ]; then
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
        "${params.align_to_ref}" \\
            "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    else
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
        -M \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
            "${params.align_to_ref}" \\
        "${fastq}" |
        samtools view \\
            -C \\
            -T params.align_to_ref \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    fi
    """
}

```

# Exact error

Here's an example of the intermittent errors I'm getting. I would have assumed that some read (e.g., read # `2_102_429` in the `.fastq`) is malformed, but as I mentioned, running `bash .command.sh` in the working directory for the failed run finishes without issue.

```
Command error:
  [M::mem_pestat] skip orientation FF
  [M::mem_process_seqs] Processed 700000 reads in 109.002 CPU sec, 16.411 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (24, 301945, 4, 10)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (118, 180, 399)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 961)
  [M::mem_pestat] mean and std.dev: (179.38, 124.36)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 1242)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (247, 319, 420)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 766)
  [M::mem_pestat] mean and std.dev: (342.90, 127.88)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 939)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (326, 764, 1452)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 3704)
  [M::mem_pestat] mean and std.dev: (893.30, 530.98)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 4830)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 97.823 CPU sec, 14.292 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (14, 196658, 9, 15)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (206, 422, 1021)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2651)
  [M::mem_pestat] mean and std.dev: (506.92, 458.91)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3466)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (248, 322, 429)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 791)
  [M::mem_pestat] mean and std.dev: (345.17, 129.22)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 972)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (485, 777, 1315)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2975)
  [M::mem_pestat] mean and std.dev: (961.60, 569.27)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3805)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 275.897 CPU sec, 39.910 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [W::sam_read1_sam] Parse error at line 2102429
  samtools view: error reading file "-"

```

A previous sample complained about an `@SQ` header line missing the `LN:` tag. It's like the pipe gets interrupted and lines get truncated, or something.

```
  [E::sam_hrecs_update_hashes] Header includes @SQ line "KI270330" with no LN: tag
  samtools view: failed to add PG line to the header

```

Anyone seen this before?

# Answer
I haven't seen that particular error before, but I'd first make sure your pipeline errors aren't being masked by your shell. Consider adding the following to your `nextflow.config`:

```
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

```

IIRC, the default shell (as provided by Nextflow) does not include the pipefail option for pipelines.

---

The only other thing I can think of is to make sure your reference FASTA (and BWA index files) are localized in the workDir. You might find the intermittent (filesystem?) errors maybe go away even if you are staging using symlinks.

```
params.ref_fasta = './bwa_index/ref.fasta'

align_to_ref = Channel
    .fromPath( [ params.ref_fasta, "${params.ref_fasta}.*" ] )
    .collect()

process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)
    path align_to_ref

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    if( params.output_format == "bam" ) {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -1 \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    } else {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -C \\
            -T "${align_to_ref.first()}" \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    }
}

```

Note that in the first case (where `params.output_format == "bam"`), the reference FASTA file itself is not required by BWA MEM. Only the BWA index files are needed here.
> 2 votes

---
Tags: sequence-alignment, samtools, bwa, nextflow
---


# Question
Title: Tajima-Nei Distance estimate with BioPerl
I have been trying to estimate Tajima-Nei distance for my data (link below). The output is s phylogenetics distance (pairwise) matrix which accommodates a detailed transition transversion and the Jukes-Cantor correction.

I´m following this protocol from BioPerl: https://metacpan.org/pod/Bio::Align::DNAStatistics

I have 314 sequences in a fasta file and another file with the list of IDs. Fasta:

```
>AVP78031.1
----atgttgtttttcttgtttcttcagttcgccttagtaaactc---------------
------------------------------------ccagtgtgttaacttgacaggcag
a----------------accccactcaatcccaattat--actaattcttcacaaagagg
...

```

IDs:

```
AVP78031.1
...

```

I´m using this Perl script to calculate the Tajima-Nei Distance in a pairwise comparison (314 * 314):

```
use strict;
use warnings;
use Bio::AlignIO;
use Bio::Align::DNAStatistics;

my $file = $ARGV[0];
my $idfile = $ARGV[1];

if ($file eq "" ) {
  $file = "NT_MSA_S_protein.fasta";
} elsif ($idfile eq "" ) {
  $idfile = "NT_ID_S_protein.csv";
}

#### Considerando un archivo
my @contentIDS;

open (LIST, $idfile) or die;
while (my $l = <LIST>) {
  $l =~ s/\n//g; # delete newline
    $l =~ s/\r//g; # delete CR
  next if (length($l) < 1);
  push @contentIDS, $l;
}
close LIST;

#### .... IDs list
my $stats = Bio::Align::DNAStatistics->new();
my $alignin = Bio::AlignIO->new(-format => 'fasta', -file   => $file);  ### $file: MSA file
while (my $aln = $alignin->next_aln) {
  #print "reading...A\n"; ### DIAG
  my $matrix = $stats->distance(-align => $aln, -method => 'Tajima-Nei');
  #print "reading...B\n"; ### DIAG
  ### Obtaining values for each pair (DISTANCE!)
  WL1:
  foreach my $aaa (@contentIDS) { ### ID #1
    WL2:
    foreach my $baa (@contentIDS) { ### ID #2
    next (WL2) if ($aaa eq $baa);
  my $data =  $matrix->get_entry($aaa, $baa);
  #($data = 0) if ($data < 0);
    print "DISTANCE\t$aaa\t$baa\t$data\n";
    } # END WL2
  } # END WL1
}

exit;

#

```

This script work it when I tried with small data, however, when I try with my real data this is the error message

```
MSG: Must provide a DNA alignment to Bio::Align::DNAStatistics, you provided a protein
---------------------------------------------------
Can't locate object method "get_entry" via package "0" (perhaps you forgot to load "0"?) at Tajima-Nei_Distance_NV.pl lin$

```

This is weird because I review my data for ambiguous characters and the characters are in majority "atcg" and on some occasions "n", at least that there are other ambiguous characters (maybe) that represent protein sequence. I really don´t understand the message because the fasta file is clearly a nucleotide sequence.

Anyone with experience using Bioperl and estimation Tajima-Nei distance?

Link: https://github.com/MauriAndresMU1313/Example\_Tajima-Nei\_Distance\_Bioperl/tree/main Any comment or help is welcome!! Thank!

# Answer
The problem is the code itself.

I suspect there are R (argenine) and Y (tyrosine) in the sequence data. These are common ambiguities but also represent amino acids. The distance estimate has not accommodated these nucleotide ambiguities but is doing a getter setter for amino acids.

It would be easier to perform the same calculation in Mega X than debugging the code.
> 1 votes

---
Tags: statistics, phylogeny, phylogenetics, perl
---