# Question
Title: BAM files with no RNAME and POS, how to map contents to SNPs?
I have a set of 4 .bam files containing the exome of an individual, around 400 MB each. I used samtools to generate a 2.4 GB .sam file out of one of the .bam files, and I found it contains lines with 101 base sequences but seemingly no RNAME or POS information on them. The first columns of the first lines in the read like this (I won't include here the SEQ, QUAL and optional fields):

```
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:18818    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10004:32033    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:12712    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:21261    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:34882    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10022:07764    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:11569    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10031:26537    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:21010    141 *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    77  *   0   0   *   *   0   0
210713_A00776_0149_AHCN5JDSX2:1:1101:10041:27085    141 *   0   0   *   *   0   0

```

Most fields are nullified. As far as I can understand, subsequent lines with the 77 and 141 flags form corresponding pairs of bases and the QNAME field first element is an individual identifier, as it won't change through the head and tail of the file. Same goes for the second element (1). The last three elements change and might hold the information I need, but I'm not sure how to read them.

**What I want to do is to map bases to SNPs (or a nucleotide identifier like NC\_000006.11).**

I have no real experience with bioinformatics besides enriching a file from MyHeritage (which contains around 609,000 base pairs with a corresponding rSNP number and chromossome) and crossed it with a reference human genome database I downloaded from NBCI. As it's a personal project, I *somewhat* validated the final result by correctly spotting a few known mutations/polymorphisms for that individual tested through usual lab tests, like the one output by the code I've wrote:

```
 RS1799945 (Chromosome 6 ): CG (NC_000006.11:g.26091179C>G, snv), https://www.ncbi.nlm.nih.gov/snp/rs1799945
          Reference Allele: C
                 Frequency:  9.6% (0.09631744626441) (~ 1:10)
              Loci (Genes): HFE  (homeostatic iron regulator)
                          : LOC108783645  (HFE antisense RNA)
        Detected Mutations: c.187C>G
             Variation URL: https://www.ncbi.nlm.nih.gov/clinvar/variation/10/
                  Diseases: not provided (conflicting-interpretations-of-pathogenicity), None (other), Hereditary hemochromatosis (pathogenic), Cystic fibrosis (risk-factor), Alzheimer disease (pathogenic), Familial porphyria cutanea tarda (None), Hemochromatosis type 1 (None), Microvascular complications of diabetes 7 (None), Transferrin serum leve, Microvascular complications of diabetes 7 (risk-factor), Hemochromatosis type 1 (pathogenic), Hemochromatosis type 2 (pathogenic)
                    PubMed: https://www.ncbi.nlm.nih.gov/pubmed/8696333, https://www.ncbi.nlm.nih.gov/pubmed/8896550, https://www.ncbi.nlm.nih.gov/pubmed/8943161, https://www.ncbi.nlm.nih.gov/pubmed/9106528, https://www.ncbi.nlm.nih.gov/pubmed/9326341, https://www.ncbi.nlm.nih.gov/pubmed/9462220, https://www.ncbi.nlm.nih.gov/pubmed/10673304, https://www.ncbi.nlm.nih.gov/pubmed/11358905, https://www.ncbi.nlm.nih.gov/pubmed/11423500, https://www.ncbi.nlm.nih.gov/pubmed/11812557, https://www.ncbi.nlm.nih.gov/pubmed/11874997, https://www.ncbi.nlm.nih.gov/pubmed/11904676, https://www.ncbi.nlm.nih.gov/pubmed/12429850, https://www.ncbi.nlm.nih.gov/pubmed/14673107, https://www.ncbi.nlm.nih.gov/pubmed/15347835, https://www.ncbi.nlm.nih.gov/pubmed/15858186, https://www.ncbi.nlm.nih.gov/pubmed/18194558, https://www.ncbi.nlm.nih.gov/pubmed/18199861, https://www.ncbi.nlm.nih.gov/pubmed/18566337, https://www.ncbi.nlm.nih.gov/pubmed/18603647, https://www.ncbi.nlm.nih.gov/pubmed/18795173, https://www.ncbi.nlm.nih.gov/pubmed/19084217, https://www.ncbi.nlm.nih.gov/pubmed/19159930, https://www.ncbi.nlm.nih.gov/pubmed/19165391, https://www.ncbi.nlm.nih.gov/pubmed/19237423, https://www.ncbi.nlm.nih.gov/pubmed/19401444, https://www.ncbi.nlm.nih.gov/pubmed/19554541, https://www.ncbi.nlm.nih.gov/pubmed/19862010, https://www.ncbi.nlm.nih.gov/pubmed/19879291, https://www.ncbi.nlm.nih.gov/pubmed/19884647, https://www.ncbi.nlm.nih.gov/pubmed/20110814, https://www.ncbi.nlm.nih.gov/pubmed/20301613, https://www.ncbi.nlm.nih.gov/pubmed/20435227, https://www.ncbi.nlm.nih.gov/pubmed/20556870, https://www.ncbi.nlm.nih.gov/pubmed/20659343, https://www.ncbi.nlm.nih.gov/pubmed/20846924, https://www.ncbi.nlm.nih.gov/pubmed/20876667, https://www.ncbi.nlm.nih.gov/pubmed/21138591, https://www.ncbi.nlm.nih.gov/pubmed/21149283, https://www.ncbi.nlm.nih.gov/pubmed/21208937, https://www.ncbi.nlm.nih.gov/pubmed/21540342, https://www.ncbi.nlm.nih.gov/pubmed/21553046, https://www.ncbi.nlm.nih.gov/pubmed/21679129, https://www.ncbi.nlm.nih.gov/pubmed/21909115, https://www.ncbi.nlm.nih.gov/pubmed/21978626, https://www.ncbi.nlm.nih.gov/pubmed/22074419, https://www.ncbi.nlm.nih.gov/pubmed/22100073, https://www.ncbi.nlm.nih.gov/pubmed/22125427, https://www.ncbi.nlm.nih.gov/pubmed/22232660, https://www.ncbi.nlm.nih.gov/pubmed/22408404, https://www.ncbi.nlm.nih.gov/pubmed/22611049, https://www.ncbi.nlm.nih.gov/pubmed/22735619, https://www.ncbi.nlm.nih.gov/pubmed/22815867, https://www.ncbi.nlm.nih.gov/pubmed/23389292, https://www.ncbi.nlm.nih.gov/pubmed/23446634, https://www.ncbi.nlm.nih.gov/pubmed/23792061, https://www.ncbi.nlm.nih.gov/pubmed/23935582, https://www.ncbi.nlm.nih.gov/pubmed/23996192, https://www.ncbi.nlm.nih.gov/pubmed/24033266, https://www.ncbi.nlm.nih.gov/pubmed/24121126, https://www.ncbi.nlm.nih.gov/pubmed/24398642, https://www.ncbi.nlm.nih.gov/pubmed/24755831, https://www.ncbi.nlm.nih.gov/pubmed/24798971, https://www.ncbi.nlm.nih.gov/pubmed/24825461, https://www.ncbi.nlm.nih.gov/pubmed/25071582, https://www.ncbi.nlm.nih.gov/pubmed/25112650, https://www.ncbi.nlm.nih.gov/pubmed/25287020, https://www.ncbi.nlm.nih.gov/pubmed/25741868, https://www.ncbi.nlm.nih.gov/pubmed/26024779, https://www.ncbi.nlm.nih.gov/pubmed/26025379, https://www.ncbi.nlm.nih.gov/pubmed/26159428, https://www.ncbi.nlm.nih.gov/pubmed/26365338, https://www.ncbi.nlm.nih.gov/pubmed/26416403, https://www.ncbi.nlm.nih.gov/pubmed/26597663, https://www.ncbi.nlm.nih.gov/pubmed/26695521, https://www.ncbi.nlm.nih.gov/pubmed/26716680, https://www.ncbi.nlm.nih.gov/pubmed/26833210, https://www.ncbi.nlm.nih.gov/pubmed/27221532, https://www.ncbi.nlm.nih.gov/pubmed/27255824, https://www.ncbi.nlm.nih.gov/pubmed/27280446, https://www.ncbi.nlm.nih.gov/pubmed/27317329, https://www.ncbi.nlm.nih.gov/pubmed/27437086, https://www.ncbi.nlm.nih.gov/pubmed/27547017, https://www.ncbi.nlm.nih.gov/pubmed/27584680, https://www.ncbi.nlm.nih.gov/pubmed/27657935, https://www.ncbi.nlm.nih.gov/pubmed/27661980, https://www.ncbi.nlm.nih.gov/pubmed/28151915, https://www.ncbi.nlm.nih.gov/pubmed/28684612, https://www.ncbi.nlm.nih.gov/pubmed/28820869, https://www.ncbi.nlm.nih.gov/pubmed/29167213, https://www.ncbi.nlm.nih.gov/pubmed/29201641, https://www.ncbi.nlm.nih.gov/pubmed/30409984, https://www.ncbi.nlm.nih.gov/pubmed/30536387, https://www.ncbi.nlm.nih.gov/pubmed/30798813, https://www.ncbi.nlm.nih.gov/pubmed/31220083, https://www.ncbi.nlm.nih.gov/pubmed/31226389, https://www.ncbi.nlm.nih.gov/pubmed/31469255, https://www.ncbi.nlm.nih.gov/pubmed/31970519
             Polymorphisms: c.118=, c.187=, c.118C>G, c.187C>G, g.8671=, g.8671C>G, g.26090951=, g.26091179=, g.26090951C>G, g.26091179C>G, p.His40=, p.His63=, p.His40Asp, p.His63Asp, n.84=, n.281=, n.309=, n.84G>C, n.281C>G, n.309C>G
     Clinical Significance: RCV000000026.18, RCV000000027.5, RCV000175607.6, RCV000394716.6, RCV000763144.1, RCV000844708.1, RCV000991133.1
                      HGVS:

```

*This question was also asked on Reddit*

# Answer
> 2 votes
These lines (with no REF or POS) are unmapped reads. If those lines are the entirety of the BAM file, then it is simply another file format to store the raw sequenced reads (as an alternative to FASTQ files). It's somewhat common for BAM files to include *all* reads, with mapped reads sorted by chromosome and position at the start, and unmapped reads appearing at the end of the file. This means that only the BAM file needs to be transferred and stored, rather than a BAM file and a separate 'unmapped' file (which can change size/content depending on the reference genome and mapper used).

FWIW, the lines of a BAM file can be viewed using `samtools view`; there's no need to convert to SAM first:

```
samtools view -h input.bam | less

```

If you only want to look at the unmapped reads, then you can filter for inclusion using the read flag 4:

```
samtools view -f 0x04 -h input.bam | less

```

If you want to *exclude* unmapped reads, you can filter that flag for exclusion (capital `F` instead of lower-case `f`):

```
samtools view -F 0x04 -h input.bam | less

```

# Answer
> 0 votes
This is an unmapped .bam file. Someone just converted the fastq to bam format, but no one has aligned these to anything. So you need to find an aligner program and a genome and align them.

---
Tags: sam, snp, samtools, identifiers
---


# Question
Title: Perform protein structure-based sequence alignment in Python
I am looking for a Python package that performs pairwise structural alignment of protein structures (i.e., PDB files) and returns a sequence alignment. PyMOL is able to do this through the GUI, for example:

For two protein PDBs, one can be aligned to the other with the "align" menu or `align` command:

With it, the sequences appear aligned in the sequence viewer according to the structural alignment, and an "alignment" object is created, which can be saved with the command:

`save A_aligned_with_B.aln, aln_structure_1_to_structure_0`

which produces the `A_aligned_with_B.aln` alignment file in Clustal format (nevermind the ?s, which correspond to special HIS resnames):

```
CLUSTAL

structure_0  YIYQDSISLPWKVLLVMLLALITLATTLSNAFVIATVYRTRKL?TPANYLIASLAVTDLLVSI
structure_1  EEQGNKL?WAALLILMVIIPTIG-----GNTLVILAVSLEKKLQYATNYFLMSLAVADLLVGL
                                         .*..**..*....**....**...****.****..

structure_0  LVMPISTMYTVT-GRWTLGQVVCDFWLSSDITCCTASIL?LCVIALDRYWAITDAVEYSAKRT
structure_1  FVMPIALLTIMFEAMWPLPLVLCPAWLFLDVLFSTASIM?LCAISVDRYIAIKKPIQANQYNS
             .****......  ..*.*..*.*..**..*....****.***.*..***.**........ ..

structure_0  PKRAAVMIALVWVFSISISLPP-FFWRQAKAEEEVSE-CVVNTD?IL-YTVYSTVGAFYFPTL
structure_1  RATAFIKITVVWLISIGIAIPVPIKGIETDVDNPNNITCVLTKERFGDFMLFGSLAAFFTPLA
             ...*...*..**..**.*..*  .              **....    ........**..*..

structure_0  LLIALYGRIYVEARSRILKQTPNRTGKLLEKKKLMAAR--ERKATKTLGIILGAFIVCWLPFF
structure_1  IMIVTYFLTI?ALQKKAYLVKNKPPQRLTIGKKSVQTISNEQRASKVLGIVFFLFLLMWCPFF
             ..*..*....                              * .*.*.***....*...*.***

structure_0  IISLVMPICKDACW--F?LAIFDFFTWLGYLNSLINPIIYTMSNEDFKQAF?KLIRFK--
structure_1  ITNITLVLCDS-CNQTTLQMLLEIFVWIGYVSSGVNPLVYTLFNKTFRDAFGRYITCNYR
             *.......*   *   .  .....*.*.**..*..**..**..*..*..**

```

I need a Python package that is easily installable/distributable that accepts two PDBs or protein structures and returns the structural alignment of the sequences somehow. All within a Python script, and not creating intermediate files, e.g. for the alignment.

There are plenty of Python tools that allow structural superimposition and RMSD calculation (e.g., Biopython), but I haven't found any that are able to return a sequence alignment from the structural alignment.

**Edit:**

Possible answer, thanks to the answer and comments from Matteo Ferla.

```
import pymol2
with pymol2.PyMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar', object='aln')
    aln = pymol2.cmd2.pymol.exporting.get_alnstr('aln', _self=pymol.cmd) # I don't know if this line could be cleaner
    # aln is itself a string

```

# Answer
> 2 votes
PyMOL in itself can be used as python package. The GUI app has its own Python, but in a conda environment you can install it as a package. Any PyMOL app command can be used as a cmd command. In PyMOL Wiki you see both. They follow a similar pattern and actually have more functionality.

```
import pymol2
with pymol2.PYMOL() as pymol:
    pymol.cmd.load('foo.pdb', 'foo')
    pymol.cmd.load('bar.pdb', 'bar')
    pymol.cmd.align('foo', 'bar')
    pymol.cmd.save('foobar.aln')

```

For more discussion of how the alignment output object is used see Assessing PyMol sequence alignment object

---
Tags: sequence-alignment, proteins, protein-structure
---


title: syndip dataset for benchmark variant

text:
# Question
Title: syndip dataset for benchmark variant
I have a question about syndip dataset : https://github.com/lh3/CHM-eval . I'm struggling to find the syndip vcf.

In the release ( https://github.com/lh3/CHM-eval/releases ), we have a file named : `rep2.37.broad.hc.raw.vcf.gz`, that i don't know what it is. And we have a file named `CHM-evalkit-20180222.tar` which contains `full.37m.vcf` and other files (bed, eval ...). So i did my search and according to this file: they mentioned that full.37m.vcf is the truth dataset. ( https://www.biorxiv.org/content/biorxiv/suppl/2018/10/29/456103.DC1/456103-1.pdf Page 16).

The problem is that the file `rep2.37.broad.hc.raw.vcf.gz` contain variants with `MQ, DP, GQ` ... that I need to extract. But the `full.37m.vcf` doesn't contain this information.. ( just `Chrom pos ref alt and QUAL`.)

So I tried to intersect `rep2.37.broad.hc.raw.vcf.gz` with `full.37m.vcf` and take the variant that present in two files, with the DP MQ GQ in `rep2.37.broad.hc.raw.vcf.gz`.

Is that okay? Since I don't know what is `rep2.37.broad.hc.raw.vcf.gz`.

And I also noticed that the QUAL in the `full.37m.vcf` is always 30. Is it normal?

# Answer
> 1 votes
If you want the truth VCF, it is `full.37m.vcf`. `rep2.37.broad.hc.raw.vcf.gz` is an example call set by GATK.

---
Tags: vcf, public-dataset, variants
---


title: recover files after groupTuple

text:
# Question
Title: recover files after groupTuple
I have a process that emits BAM files with a key containing the sample name:

```
process source{
  output:
    tuple val(sample), path("sample.bam") into ch

  [...]

}

```

1st question: I have several files per `sample`, let's assume 2 for now. I want to merge them together with `samtools merge`. I can regroup that channel with `groupTuple()`, but how can I recover them in the next process? Considering they all have the same name:

```
grouped_ch = ch.groupTuple(by: 0, size: 2)

process merge{
  input:
    grouped_ch

  shell:
  '''
  samtool merge sample.bam sample.bam > !{sample}.bam
  '''
}

```

Is there perhaps a way in the `input:` block to say that I expect a list, and name the list components?

2nd question: In practice, I have a variable number of "sample.bam" files for each `sample` value, does that change the way to define the `input:`? It also appears I will need to use the built-in `groupKey()`, is there documentation available? I have read this Github issue, but it would help if a general description was available.

# Answer
My preference is to avoid working with files with the same name as much as possible. However, sometimes this might be unavoidable. Fortunately, file name collisions can be avoided when working with multiple input files by using the `*` and `?` wildcards:

> When a target file name is defined in the input parameter and a collection of files is received by the process, the file name will be appended by a numerical suffix representing its ordinal position in the list.
>
> The target input file name can contain the * and ? wildcards, that can be used to control the name of staged files.

Note that there is also the `stageAs` path option (see input of type 'path'), which can also take a name pattern (using the above wildcards) to let you use a variable in your script. For example:

```
process samtools_merge {

    conda 'samtools=1.15.1'
    debug true
    stageInMode 'rellink'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    ls -nl
    """
}

workflow {

    Channel.fromPath( ['./path/to/files/*.bam', './path/to/more_files/*.bam'] ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple(by: 0, size: 2) \
        | samtools_merge
}

```

Results:

```
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [elated_montalcini] DSL2 - revision: 1427d64a1d
executor >  local (3)
[83/25654c] process > samtools_merge (3) [100%] 3 of 3 ✔
total 184
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample1.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample1.bam
-rw-r--r-- 1 1000 985 187911 May 25 12:18 sample1.bam

total 128
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample2.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample2.bam
-rw-r--r-- 1 1000 985 129947 May 25 12:18 sample2.bam

total 196
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample3.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample3.bam
-rw-r--r-- 1 1000 985 197225 May 25 12:18 sample3.bam

```

If each group contains a variable number of files, you can keep the same input declaration (as above) but you will want use a special groupKey object so that the collected values can be streamed as soon as possible. It's not well documented, but here's a simple example that sorts the input BAM files prior to merging. Once all of the 'samtools\_sort' jobs are finished for a particular sample, that sample can proceed to merging without us having to wait for all the other 'samtools\_sort' jobs to be finished (which would otherwise be the case if the 'size' parameter had not been specified):

```
process samtools_sort {

    tag { bam.baseName }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam)

    output:
    tuple val(sample), path("${bam.baseName}.sorted.bam")

    """
    samtools sort -o "${bam.baseName}.sorted.bam" "${bam}"
    """
}

process samtools_merge {

    tag { sample }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    """
}

```

```
workflow {

    Channel.fromPath( params.bam_files ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple() \
        | map { sample, files -> tuple( groupKey(sample, files.size()), files ) } \
        | transpose() \
        | set { bam_files }

    samtools_sort( bam_files ) \
        | groupTuple() \
        | samtools_merge \
        | view()
}

```

Results:

```
$ nextflow run script.nf -ansi-log false --bam_files './path/to/*/*.bam'
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [cheesy_swanson] DSL2 - revision: da1a82a059
[04/643129] Submitted process > samtools_sort (sample3)
[1a/93a939] Submitted process > samtools_sort (sample1)
[93/fde0a0] Submitted process > samtools_sort (sample3)
[db/208bac] Submitted process > samtools_sort (sample3)
[8f/f2a231] Submitted process > samtools_sort (sample2)
[f3/8f0f9f] Submitted process > samtools_sort (sample2)
[bb/cca23e] Submitted process > samtools_merge (sample1)
[sample1, /home/steve/testing/work/bb/cca23e016a641b201b7e42f9db74ef/sample1.bam]
[ff/bd3a2a] Submitted process > samtools_merge (sample2)
[sample2, /home/steve/testing/work/ff/bd3a2ac83a0691d2b80b04ef97b755/sample2.bam]
[94/a4f8eb] Submitted process > samtools_merge (sample3)
[sample3, /home/steve/testing/work/94/a4f8eb483602e9335b28f5d40058c8/sample3.bam]

```

Note that the val specified as input to 'samtools\_sort' is an object of class `nextflow.extension.GroupKey`. In the example above, the input to 'samtools\_merge' is also a groupKey - but it doesn't need to be. If you need a string value for whatever reason, you can just call '.toString()' to get it back:

```
    ...

    samtools_sort( bam_files ) \
        | groupTuple() \
        | map { sample_key, files -> tuple( sample_key.toString(), files ) } \
        | samtools_merge \
        | view()

```
> 3 votes

---
Tags: nextflow
---


# Question
Title: Looking for a good tool for repeat masking in long sequences
The repeat masker site has a size limit of 90kb. Where can I mask repeats in a 500kb sequence?

Ideally, I need a thorough repeat masking of known and de-novo found repeats of various complexity. I am working on Windows but can access Unix too.

# Answer
> 4 votes
RepeatMasker can downloaded an run locally. There's then no size limit, though do note that it's quite slow for mammalian-sized sequences (you'll want to split it into chunks and run in parallel in such cases). For a small 500kb sequence that won't pose an issue.

# Answer
> 3 votes
RepeatFinder may suffice by Volfovsky et al. (2001) and is available for download.

# Answer
> 2 votes
From: http://meme-suite.org/doc/glam2\_tut.html, here's a few tools that do various repeat masking tasks.

> These tools mask various kinds of repeats: Repeatmasker, Censor, seg/pseg/nseg, DUST, XNU, DSR, GBA.

Censor, seg/pseg/nseg, dust, DSR all have locally available applications. DSR, GBA, and CAST are algorithms for repeat masking.

# Answer
> 0 votes
Please use this python program to mask elements -

```
dna = input("Please enter the DNA sequence or chromosome whose repetitive elements need to be masked:")
repele = ["aaaa", "tttt", "cccc", "gggg"]
counter=0
str1=str(dna)
for i in repele:
    str1 = str1.replace(repele[counter], "x")
    counter=counter+1
print(str1)

```

---
Tags: repeat-elements, repeatmasker, repeat
---


# Question
Title: Generate consensus protein sequence from relatively gappy alignment?
I would like to know what would be the best way to generate a full-length consensus sequence?

I want to obtain a single representative sequence from a relatively gappy multiple sequence alignment of peptides.

See example below.

So far I have attempted `Consensus.py` from https://github.com/jrjhealey/bioinfo-tools but it just gives me gaps, even with a threshold close to zero.

```
>HX5_33766:7:1222:1022:72139
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1211:6319:29841
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1222:1022:72139.2
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1219:12723:20928.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1119:11069:37612
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1219:12723:20928
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1211:6319:29841.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1219:12723:20928.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------T-ISTSYKDIQVNGGELILIHQNPGGFCAY--------------
>HX5_33766:7:1107:2564:39757
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1211:6319:29841.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIIINCDGY------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1119:11069:37612.2
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165.2
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1106:31832:24286
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:10013:25165.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMK-------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314.2
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1212:30371:5686.1
----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------CENEKGCIITSPSKDETMKVQDNSIXLGCN--------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1221:9364:63314.1
-------------------------------GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1119:11069:37612.1
----------------------------------------------------------------------------------------------------------------------------------------VLLGAAIIQGLGLLLC--F----------------------TYI----------------------CLHFYASQ-V----------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1107:2564:39757.1
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1106:31832:24286.1
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881.1
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1222:1022:72139.1
---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------SEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFKDKVYFNVT--------------------------------------------
>HX5_33766:7:1107:2564:39757.2
-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------VFSLKGYFSEEL---SLSLYYRKGRGPLF-------SLSKVTSVDSIGVAYLAFK----------------------------------------------------
>HX5_33766:7:1106:31832:24286.2
----------------------------------------------------------------------SVYLRL--IVFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1117:3691:48881.2
-------------------------------------------------------------------------------VFP-PSPPLWRRG---QSP------SGLLPFATV-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213.1
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1121:4320:48213.2
MTVISQSLKSELECHKNNSLLAQ-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1214:16701:50287
----------------------------------------------------------------------------------------------------------------------------------------------R--KQR--EXPKIEQLSAFTRSRSPESAFSNNHSVTYFHTLWGD------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
>HX5_33766:7:1220:24545:3085
------------------------------------------------------------------------------------------------------------------------------------------------------------------RSRSPESAFSNNHSVTYFHTLWGD-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------R

```

# Answer
Given the alignment you've presented, there are obvious groups of sequence similarity, but I'm seeing a lot of different peptide sequences in this (e.g. `CENEKGCIITSPSKDETMKVQDNSIXLGCN` and `GKFRLVTKKFPRLPAEAAVLKLRECVLLFQLPLISAPGTSVYLR` as two obviously different peptides). Have you done a pre-filter for similarity to make sure you're looking at the same proteins for alignment?

Multiple alignment will not work properly with peptides that are very different from each other. You need to group proteins based on similarity first (e.g. with Linclust), then carry out multiple alignments within the groups.
> 1 votes

# Answer
```
from Bio import AlignIO
from pathlib import Path
from Bio.Align import AlignInfo

path = '/my/path/'
test1 = 'testset.fa'
mypath = Path(path,test1)

myalign = AlignIO.read(mypath, "fasta")
summary = AlignInfo.SummaryInfo(myalign)
consensus = summary.dumb_consensus()
print (consensus)

```

This will give a single consensus and 'X' low prevalence ambiguity. Easy to check whether its what is required because the code is ready to go.
> 0 votes

---
Tags: fasta, phylogenetics, proteins, multiple-sequence-alignment
---


# Question
Title: Phylogenetic tree building from proGenomes database for shotgun metagenomics
For some weeks I'm fighting with an issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account phylogenetic distance.

I’m trying to create a phylogenetic tree from the representative genomes of proGenomes v.2.1 database. However, I do not have a representative genome for each taxa id. In contrast, I have a multifasta files for more than 100 contigs for each taxa ID (file freeze12.contigs.representatives.fasta.gz).

For each taxaID (the first numbers in each header of the multifasta file is NCBI taxon ID), I would like to bin the contigs into the biggest scaffold possible, and then perform an alignment (through MUSCLE as examples) followed by RaXML-NG to construct a Newick tree.

However, I cannot find any software for binning contigs with only the contigs as inputs. MetaBAT needs bam files from shotgun read mappings to the contigs in question (not present in the proGenomes database). I found another software (Scaffold\_builder), but it needs a representative reference genome. And in microbiology I have more than 100 complete genomes for each NCBI tax id...

For this reason, I ask here an open question if exists some software that can perform of contigs taking into account only a multifasta file of contigs?

# Answer
I would not recommend the approach that you describe. I am not sure that I completely understand it. For example, metagenomic binning by e.g. MetaBAT is not a scaffolding operation, which returns a single sequence. Rather, they are methods to separate a pool of contigs into groups with similar sequence characteristics indicating they come from the same genome.

It sounds like you already have a group of contigs from each genome, you just want them to all be part of a single sequence in a FASTA file, as would be suggested by your interest in the reference-guided scaffolding approach scaffold\_builder.

I would propose that you instead use a tool like GTDBtk to place each of your genomes (each represented by a multifasta of contigs, which is the expected input) into an existing taxonomic tree, and then use their placements in that tree as a model of the taxonomy of your genomes, rather than trying to infer it de novo. I assure you that it is a high quality phylogeny that no one will criticize, and it avoids the possibility of introducing errors into your phylogeny through problematic scaffolding operations.

I would suggest reading GTDBtk documentation for more information.

## A narrower answer to the original question

I realized that there is an answer to the question of binning contigs without any additional information. I do not recommend it as I think it would give you positively wrong answers because you are starting with your genomes already binned, and it would not provide you with scaffolds, which I think is what you actually want. But there is a way to perform this task.

You can compute tetranucleotide frequencies from contigs and use them to cluster arbitrary sets of contigs into genomes or other groupings based on similarities. This has a significant rate of error, but it does more or less work. It's an older approach.
> 1 votes

# Answer
Thank you another time @Maximilian Press. Yes, GTDB-Tk it's a robust software but although the merging of archaea and bacteria taxa is still a current research topic, I have succeed with another tool that permits me obtaining an unique tree (including archea and bacteria).

If I work with GTDB-Tk I've obtained one tree for archaea and other for bacteria. Then to use it in phyloseq I must work with two independent phyloseq objects which is nice also but my purpose of building a tree from shotgun is only to permits me to calculate diversities based on phylogeny.

As alternative that works it was the following: to construct a phylogenetic tree from Kaiju I used the 40 marker proteins that is based on the proGenomes database used by kaiju (https://progenomes.embl.de/data/proGenomes2.1\_markerGenes.tar.gz) and then I've used the ETE-Toolkit(http://etetoolkit.org/) software to construct the tree in Newick format that I can integrate into a phyloseq object.

I post here if someone can help my strategy. It has been no easy (unless for a dummie as me) to construct a phylogenetic tree from Kaiju's taxonomic classifier.
> 0 votes

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


# Question
Title: Is this RNA seq data good based on the MDS plot
I am analyzing some RNA sequencing data from a collaborator where the effect of some ligands on the transcriptome is being looked at. I am using DESeq2 for my analysis. I am looking to compare each of the ligand treatments vs the DMSO. I have attached the MDS plot with this post. There seems to be a lot of variation between the 3 replicates across all the treatments. Based on this plot, is it worth it to use this data to run a differential expression?

It is also at this link below. https://i.stack.imgur.com/ZXhk5.jpg

# Answer
> 2 votes
You definitely see separation by treatment in PC1 and PC2, so that is an indication that the treatments are doing something. DMSO and TCDD cluster close, so these are probably similar. Still, you see a variation along PC1 and it seems that there is within-group separation in all groups that could indicate a batch effect. Check whether the samples per group that are more leftish have been processed as one batch, if so one could include that information into the design, see DESeq2 vignette on batch effects. Alternatively, check tools such as RUVseq or SVA to detect unknown batch effects to correct for this within-grup separation.

---
Tags: rna-seq, pca
---


# Question
Title: How can I improve or otherwise investigate an unreliable genome tree?
**Summary** My genome tree doesn't agree with my gene trees and I get the feeling that my genome tree might be wrong, possibly due to long branch attraction, but I don't know how to check/fix it.

**Background** I have a set of genes of interest recovered from genomes and metagenomes. I also have some bins built from the metagenomes of interest.

**Previous work** I've built multiple gene trees, which tend to agree with each other, as well as with other phylogenetic methods (16S rRNA, initial genome trees where available, conserved orientation of the genes of interest). So at this point I though I had a pretty good idea of how my genes / species of interest are linked to the species tree.

**Concern** I've assembled the metagenomes and recovered bins. I've refined bins of interest with anvio, based on sequence composition / differential coverage, and kept only the bins which appeared quite high quality (\>70% completion\*, \<10 - often \<5% - redundancy). Thus contamination in the metagenome assemblies (MAGs) is unlikely. The contamination would have to happen multiple times in the same way to different samples from different environments, at the same coverage, to get the same effect.

* with the exception of a bin of particular interest which I accepted with 62%

I have de-replicated these bins, added some known genomes, and tried to build a species tree with Orthofinder (with the proteomes). There is a problem occurs with a clade of 3 bins I will call 'X'.

**Trees** Based on the gene-of-interest phylogenies, X looks like probably a sister-group to the rest of known genomes from order 'A'. Which is very interesting. With the 16S rRNA, it's not quite as clear (there are no 16S rRNA sequences in the bins); however, some 16S rRNA amplified from one of the original samples fit the bill as 'sister-clade to basically everything else in group A'). The first (100 or so) BLAST hits for this 16S sequence are all from group 'A', still.

**Tree conflict** However, in the proteome tree, X as well as Y (Y being a sequenced genome that I suspected was related to X based on 16S data) end up in order 'B' instead of 'A'. GTDBtk also labels X as belonging to clade 'B', not 'A'.

Now,I know that genes and genomes don't always have the same evolutionary trajectories. But multiple sets of data lead me to a different conclusion than the proteome tree, so I have to at least question whether the proteome tree is correct; what's likely to be the problem (I suspect long-branch attraction), and, if it is long-branch attraction (or anything else) how can I uncover the correct signal?

Thank you for your time!

# Answer
The difficulty in answering the question is:

* The absence of "Y" in the "Gene of interest" tree
* The sampling bias in the number of taxa mostly between the "genome tree" and the "gene of interest" and "16S" tree
* The "genome" tree is unrooted so its difficult to compare it with the "16S" tree, thus its possible that if "16S" was unrooted it falls within monophyly "B"
* The trees don't have any branch or statistical support
* We don't really know the tree building method (I haven't looked it up)
* The bacteria might have significant recombination (we don't know the bacteria).

The answer is there is a lot of room for different tree artefacts here, however it could be a biologically real result. I agree "the tree" is simply a tool which uses a series of assumptions and obtaining confidence in a tree requires a lot of analysis. Vast amounts of computational time are dedicated to obtaining a robust tree.

The principle step is to compare like for like:

* same number of taxa and
* same rooting position (or unrooted).
* second step is to obtain statistical robustness for the observation.

There is a formal process for identifying conflicting trees and it needs code - I have written stuff that does this and maybe I should release it. Essentially, you would data chunk gene by gene to produce a genome of gene trees and use e.g. ETE3 or similar to plot the tree conflicts. It ain't rocket science. Thus the amount of user time is trivial. The idea is that once the "problem trees" are identified, e.g. rapid evolution.

The one observation I would say is that the "Y" branch looks a little bit weird and this could indicate part of the data is has resulted from different evolutionary rate and that is the basis for long-branch attraction.

---

**Scenario 1** It could be that the 16S and "genome" tree are part of the same tree (monophyly B), they are just rooted differently. Thus the conflict is just a small number of genes and this could lead to ...

**Scenario 2** Keep in mind bacteria do undergo 'horizontal transmission' and if 'gene of interest' is under strong selective pressure recombination is entirely possible.
> 2 votes

# Answer
It seems like you've got the right approach, in trying out a few different genes and looking for consistency. As you say, different genes can have different selective pressures, which influences how much they change in relation to the whole genome background.

Certainly don't put too much trust in a single gene, and especially not 16S: it has a low mutation rate, can jump around and be repeated in the genome, and is an occasional target for horizontal gene transfer.

The answer may be that there *isn't* any single correct phylogeny. A primary assumption of phylogenetics is that evolutionary progression has a tree-like structure (with roots and tips), but that model is broken by horizontal gene transfer and plasmid sharing. It may be the case that you need to break the genome up into fragments, each with its own inherited history.
> 1 votes

---
Tags: phylogenetics, phylogeny, assembly, metagenome, bacteria
---


# Question
Title: How to extract residues from specific chain of interest from multiple pdb files?
I'm new to biopython & I'm trying to extract residues from specific chain of each file from a list of 316 pdb files. Eg: 1d8t has four chains A, B, C & D of which my interest is to extract just the residues of C chain. As the chain of interest differs for different files, I'm not knowing how to extract specific chain residue from each file and save it separately.

```
from Bio import *
from Bio.PDB.PDBParser import PDBParser
from Bio.PDB.Polypeptide import PPBuilder
structure = PDBParser().get_structure('1d8t', '1d8t.pdb')

model = structure[0]
chain = model['C']

for i in chain.get_residues():
    print(i)

```

With the above code I was able to extract the residues of C chains not sure how to proceed for multiple files in python.

# Answer
Create a function and run it for each pdb file.

```
from Bio.PDB import PDBParser

def print_residues(name, pdb_file, chain_key, out_file):
    parser = PDBParser()
    structure = parser.get_structure(name, pdb_file)
    model = structure[0]
    chain = model[chain_key]
    with open(out_file, 'w') as fp:
        for i in chain.get_residues():
            print(i, file=fp)

```

Then run the function for all pdb files.

```
print_residues('1d8t', '1d8t.pdb', 'C', '1d8t.txt')
# ...

```
> 1 votes

---
Tags: python, pdb, parsing
---


# Question
Title: Can Bioinformatics be used in developing a marketable product?
As we all know computer science, electrical engineering, etc. knowhows can be used in developing software and electronics/electrical products that can be sold to end-users (B2C).

Can bioinformatics be used to develop such products other than drug discovery?

Or, is it just about scientific/academic research meant for B2B clients?

# Answer
It is pivotal to vaccine research at multiple levels:

1. Deep theoretical vaccine design, which takes years to get to market (my stuff's gone into that and still stuck in advanced clinical trials).
2. Looking for variants that might have evolved immunoevasion against current vaccines involving SARS-CoV-2 (why do we need so many boosters?) and influenza.

Bioinformatic patents are tricky.

The reality is both points are not systematically automated in 'point and click' user friendliness. Its definitely of commercial value at an industrial/institutional level. However most work is submitted to the public domain and code freely available or held under a commercial firewall.

Generally, I think you underestimate how much bioinformatic work takes place in modern clinical practice and when fighting a virus how it just keeps evolving around all the obstacles thrown at it.

More generally, a lot of bioinformatics is used to direct wet-lab research because bioinformatics gets to the answers quicker.
> 1 votes

---
Tags: untagged
---


# Question
Title: Are the GIAB reference genomes (HG00N) fully phased?
Specifically, the NISTv4.2.1 HG002 dataset here (more info on Genome in a Bottle). According to their README, the HG002 v4.2.1 VCF should be fully phased.

> For HG001, HG002, and HG005, we have phased the v4.2.1 benchmark vcfs by transferring phasing information from a phased assembly-based trio-hifiasm dipcall vcf to the v4.2.1 vcf

I ask because there appear to be FAR more heterozygous variants on the second haplotype than the first (`0|1` \>\> `1|0`).

For example:

```
tail -n 500000 HG002_GRCh38_1_22_v4.2.1_benchmark.vcf | \
        cut -f 10 | cut -f1 -d: | sort | uniq -c
 302960 0/1
   1104 1/0
 190042 1/1
   2037 1/2
   3857 2/1

```

Is this simply an artifact of how these truth sets were produced? For example, if the GRCh38 (haploid) FASTA was meant to be haplotype 1, then it makes sense that it's nearly identical to HG002's haplotype 1, and haplotype 2 has more differences. Is that the case, or are only certain variants phased (such as those that pass a minimum confidence threshold, or are in a high-quality BED region)?

# Answer
> 2 votes
As the two alleles in a genotype are separated by `/` but not `|`, the main v4.2.1 VCF is not phased. The README says the phased VCFs are available in the SupplementaryFiles directory. It contains a few versions of VCFs with v4.2.1 genotype calls but with phasing from additional tools such as whatshap and hifiasm. I don't know which version is the most accurate, though.

---
Tags: vcf, phasing
---


# Question
Title: BPGA error: there was no protein data in file
I'm trying to do pan-genome analysis using BPGA and I face this error:

This error didn't show up in the other data set where all the .gbk is downloaded from NCBI. This .gbk file is our own data and has not been published in NCBI yet. There's of course protein data in it, there's this CDS qualifier "translation", but somehow BPGA didn't detect this.

The difference between our data and NCBI data is that most of the NCBI .gbk files are annotated with PGAP. On the other hand, our .gbk is annotated with Prokka (have tried it in DFAST and KBase). NCBI .gbk files that are annotated with PGAP have this CDS qualifier "protein\_id" in them. Besides that, I don't think there's any difference.

Does the BPGA script counts the protein number based on the "protein\_id" qualifier? How should I solve this? I can't annotate the .fasta file with PGAP because it's so resource-intensive. I believe that, if Prokka is installed locally, there's this "--compliant" option that can integrate the "protein\_id" qualifier into the output .gbk file. I have tried to install Prokka locally, but it didn't run well too, I encountered another error for that and all the solutions didn't work for me.

Any opinion on this issue? Thanks for reading my question.

EDIT: I tried to insert one dummy "protein\_id" on my Prokka-annotated .gbk files and it run perfectly without errors.

But when I tried to insert dummy protein\_id on my DFAST-annotated .gbk files, it won't run and faced another error. What is this "protein\_id" really defined by the script?

# Answer
@raysteven,

If you follow the documnetation of BPGA at https://sourceforge.net/p/bpgatool/wiki/Home/, you can find the example of **SUPPORTED INPUTS** for all types of input files including gbk or gbff. Unfortunately, Prokka generated files are not supported due to different or missing tags than NCBI.

One simple solution you may try is to edit your genBank file to include the tags such as **/protein\_id="Any\_unique\_protein\_id\_inside\_quotes"** and everything else shown in the gbk example.

I hope this helps.
> 1 votes

---
Tags: perl, sequence-annotation
---


# Question
Title: To create a new column with Pandas on the basis of two types of column values
I have a data frame `df`. The below table only represent top 8 rows of the data frame. There are multiple combination of `LIBRARY_NAME` and `ANCHOR_NAME`

```
BARCODE ANCHOR_NAME     LIBRARY_NAME    ANCHOR_CONC Synergy HSA     Bliss
14482   5-Fluorouracil  Cisplatin       2.5         FALSE   FALSE   FALSE
                                        10          FALSE   FALSE   FALSE
                        Docetaxel       2.5         FALSE   TRUE    TRUE
                                        10          FALSE   FALSE   FALSE
        AZD7762         Cisplatin       0.0625      TRUE    TRUE    TRUE
                                        0.25        FALSE   TRUE    TRUE
                        Docetaxel       0.0625      FALSE   FALSE   FALSE
                                        0.25        FALSE   TRUE    FALSE

```

I want to create a new column 'new\_synergy 'with pandas which is based on column `ANCHOR_CONC` and `Synergy`. There are two ANCHOR\_CONC for each anchor\_name and library\_name combination. If for any concentration of ANCHOR\_CONC the synergy is `TRUE`, for other concentration also synergy becomes `TRUE`

Eg: `new_synergy`: Condition for TRUE (for AZD7762 and Cisplatin), IF at 0.0625 `Synergy` is TRUE, it will be TRUE for 0.25 also.

`Expected data_frame`

`new_df`

```
BARCODE ANCHOR_NAME     LIBRARY_NAME    ANCHOR_CONC Synergy HSA     Bliss  new_synergy
14482   5-Fluorouracil  Cisplatin       2.5         FALSE   FALSE   FALSE    FALSE
                                        10          FALSE   FALSE   FALSE    FALSE
                        Docetaxel       2.5         FALSE   TRUE    TRUE     FALSE
                                        10          FALSE   FALSE   FALSE    FALSE
        AZD7762         Cisplatin       0.0625      TRUE    TRUE    TRUE     TRUE
                                        0.25        FALSE   TRUE    TRUE     TRUE
                        Docetaxel       0.0625      FALSE   FALSE   FALSE    FALSE
                                        0.25        FALSE   TRUE    FALSE    FALSE

```

I tried the following code

```
import pandas as pd
import numpy as np
dd = data[["BARCODE", "LIBRARY_NAME", "ANCHOR_NAME", "ANCHOR_CONC", "Synergy", "Bliss", "HSA"]]
df = pd.pivot_table(data=dd,index=["BARCODE", "ANCHOR_NAME", "LIBRARY_NAME", "ANCHOR_CONC"])
## setting conditions as pivot table step converted true to 1 and false to 0.
conditions = [(dd['Synergy'] == 1.0), (dd['Synergy'] == 0.0)]
values = ['True', 'False']
dd['new_synergy'] = np.select(conditions, values)
## However it is giving me the same results in both new_synergy and synergy column. I have no idea how to fit `ANCHOR_CONC` condition in this step.

```

# Answer
You can use groupby() as following :

`dd["new_synergy"] = dd[["ANCHOR_NAME", "LIBRARY_NAME", "Synergy"]].groupby(by=["ANCHOR_NAME", "LIBRARY_NAME"]).any()`

Explenation : you want to group your dataframe by both `ANCHOR_NAME` and `LIBRARY_NAME`, and return True if any of the `synergy`value is true, which is what `any()` does.

Note : you can check for all methods available for `GroupBy` objects here.
> 3 votes

# Answer
To elaborate on this answer, here is a full working example where I assume that `"BARCODE"`, `"ANCHOR_NAME"` and `"LIBRARY_NAME"` are used as `MultiIndex`:

Creating the test data:

```
from io import StringIO
import pandas as pd
data_str = """BARCODE\tANCHOR_NAME\tLIBRARY_NAME\tANCHOR_CONC\tSynergy\tHSA\tBliss
14482\t5-Fluorouracil\tCisplatin\t2.5\tFALSE\tFALSE\tFALSE
14482\t5-Fluorouracil\tCisplatin\t10\tFALSE\tFALSE\tFALSE
14482\t5-Fluorouracil\tDocetaxel\t2.5\tFALSE\tTRUE\tTRUE
14482\t5-Fluorouracil\tDocetaxel\t10\tFALSE\tFALSE\tFALSE
14482\tAZD7762\tCisplatin\t0.0625\tTRUE\tTRUE\tTRUE
14482\tAZD7762\tCisplatin\t0.25\tFALSE\tTRUE\tTRUE
14482\tAZD7762\tDocetaxel\t0.0625\tFALSE\tFALSE\tFALSE
14482\tAZD7762\tDocetaxel\t0.25\tFALSE\tTRUE\tFALSE
"""
# StringIO is to simulate an import from a file
data = pd.read_table(StringIO(data_str), sep="\t", index_col=[0, 1, 2])
# reset_index is here just for displaying purposes
# (MultiIndex does not seem well supported in (pandas?) markdown export)
print(data.reset_index().to_markdown(index=False))

```

| BARCODE | ANCHOR\_NAME | LIBRARY\_NAME | ANCHOR\_CONC | Synergy | HSA | Bliss |
| --- | --- | --- | --- | --- | --- | --- |
| 14482 | 5-Fluorouracil | Cisplatin | 2.5 | False | False | False |
| 14482 | 5-Fluorouracil | Cisplatin | 10 | False | False | False |
| 14482 | 5-Fluorouracil | Docetaxel | 2.5 | False | True | True |
| 14482 | 5-Fluorouracil | Docetaxel | 10 | False | False | False |
| 14482 | AZD7762 | Cisplatin | 0.0625 | True | True | True |
| 14482 | AZD7762 | Cisplatin | 0.25 | False | True | True |
| 14482 | AZD7762 | Docetaxel | 0.0625 | False | False | False |
| 14482 | AZD7762 | Docetaxel | 0.25 | False | True | False |

Creating a `DataFrame` with the new column:

```
new_synergy = data.reset_index().groupby(
    ["BARCODE", "ANCHOR_NAME", "LIBRARY_NAME"])[["Synergy"]].any()
# Rename the resulting "Synergy" column
new_synergy.columns = ["new_synergy"]
print(new_synergy.reset_index().to_markdown(index=False))

```

| BARCODE | ANCHOR\_NAME | LIBRARY\_NAME | new\_synergy |
| --- | --- | --- | --- |
| 14482 | 5-Fluorouracil | Cisplatin | False |
| 14482 | 5-Fluorouracil | Docetaxel | False |
| 14482 | AZD7762 | Cisplatin | True |
| 14482 | AZD7762 | Docetaxel | False |

Merging the two `DataFrame`s based on their indices:

```
new_data = data.merge(new_synergy, left_index=True, right_index=True)
print(new_data.reset_index().to_markdown(index=False))

```

| BARCODE | ANCHOR\_NAME | LIBRARY\_NAME | ANCHOR\_CONC | Synergy | HSA | Bliss | new\_synergy |
| --- | --- | --- | --- | --- | --- | --- | --- |
| 14482 | 5-Fluorouracil | Cisplatin | 2.5 | False | False | False | False |
| 14482 | 5-Fluorouracil | Cisplatin | 10 | False | False | False | False |
| 14482 | 5-Fluorouracil | Docetaxel | 2.5 | False | True | True | False |
| 14482 | 5-Fluorouracil | Docetaxel | 10 | False | False | False | False |
| 14482 | AZD7762 | Cisplatin | 0.0625 | True | True | True | True |
| 14482 | AZD7762 | Cisplatin | 0.25 | False | True | True | True |
| 14482 | AZD7762 | Docetaxel | 0.0625 | False | False | False | False |
| 14482 | AZD7762 | Docetaxel | 0.25 | False | True | False | False |
> 3 votes

---
Tags: python, pandas
---


# Question
Title: How to get strain names/ids contained in a multi FASTA file using seqkit?
FASTA files can be very big and unwieldy, especially if lines are at most 80 characters, one can't speed up browsing them by using `less` with `-S` to have one sequence every two lines.

How can I extract just the strain names (or sequence names, i.e. the string on the line starting with `>`) into a list? I'd like to use `seqkit` since I'm already familiar with it and know some of the flags and idioms it uses.

I googled the question but couldn't find an answer using seqkit - the examples all use `awk`, `grep`, `sed` or similar, which are alright for copy/pasting but not very memorable.

# Answer
> 6 votes
The nice thing about `seqkit` is that the answer is really simple and memorable:

```
seqkit seq -n sequences.fasta

```

This will output a text file containing one sequence name per line:

```
MPX/PT0001/2022
MPX/PT0002/2022
MPX/PT0003/2022

```

# Answer
> 4 votes
You can trivially do this with standard \*nix tools. Here are a few options:

```
sed -n 's/^>//p' file.fa
awk 'sub(/^>/,"")' file.fa
grep '^>' file.fa | tr -d '>'
perl -ne 'print if s/^>//' file.fa

```

The advantage of these is that they will work on any POSIX-compliant system and most other \*nix flavors, even if they're not POSIX such as Linux, out of the box without needing to install any external tools.

Another thing to take into account is how fast a command is, especially when dealing with large files. I ran a few tests using the solutions presented here on the human hg38 genome assembly which can be downloaded here and `seqkit` version `2.2.0` which was downloaded from here. The file was stored locally on my laptop's SSD drive, and the laptop itself is:

```
$ inxi
CPU: quad core Intel Core i7-6820HQ (-MT MCP-)
speed/min/max: 2896/800/3600 MHz Kernel: 5.17.5-arch1-1 x86_64
Up: 25d 10h 27m Mem: 15035.9/31958.1 MiB (47.0%)
Storage: 476.94 GiB (103.1% used) Procs: 452 Shell: Bash inxi: 3.3.15

```

To get more accurate results, I used a little function I have which will run the command you give it a specified number of times and then report the average *real* time across all runs. Here's the code if you want to try it:

```
timethis(){
    max=$1;
    shift;

    for com
    do
        printf '\nCOMMAND: %s\n' "$com"
    c=0;
    while [[ $c -lt $max ]]
    do
      let c++;
      ( time -p eval "$com" ) 2>&1 | grep -oP 'real.*?\K[\d\.]+'
        done |
        awk -vm=$max '{k+=$1}END{print (k/m)}';
    done
}

```

And the results. These are the average "real" (i.e. wall clock) times taken across 10 runs:

```
$ timethis 10 \
  "rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null" \
  "sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "awk 'sub(/^>/,\"\")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null" \
  "perl -nle 'push @res, $_  while s/^>//; END {print  join \"\n\",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna"

COMMAND: rg '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.474

COMMAND: grep '^>' GCA_000001405.29_GRCh38.p14_genomic.fna | tr -d '>' >/dev/null
0.583

COMMAND: sed -n 's/^>//p' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
5.62

COMMAND: perl -ne 'print if s/^>//' GCA_000001405.29_GRCh38.p14_genomic.fna  >/dev/null
6.408

COMMAND: ./seqkit seq -n GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
4.326

COMMAND: awk 'sub(/^>/,"")' GCA_000001405.29_GRCh38.p14_genomic.fna >/dev/null
13.518

COMMAND: perl -nle 'push @res, 3  while s/^>//; END {print  join "\n",@res}' GCA_000001405.29_GRCh38.p14_genomic.fna
8.455

```

As you can see above, the `grep` or `rg` and `tr` approaches only took around half a second on average, so that is clearly much faster than any of the other choices. However, `seqkit` is pretty good since it only took around 4 seconds, faster than any of the other ad hoc solutions.

# Answer
> 1 votes
Or a bit more verbose ...

```
perl -nle 'push @res, $_  while s/^>//; END {print  join "\n",@res}' myfile.fa

```

---
Tags: fasta, identifiers, seqkit, multi-fasta
---


# Question
Title: Checking that two multiline FASTAs are identical (allowing for different order and lower/upper case)
Is there a good CLI tool that can diff two multiline fasta files, ignoring order, and comparing sequences with identical names with each other pairwise.

Bonus points if lower/upper case is tolerated (as in, it's considered identical) and if also indels can be tolerated, maybe through a flag?

I could script this in Python myself, but I was wondering whether there was a custom tool for this purpose already out there. Maybe something more performant than a Biopython script.

If there's no custom tool, I'd also be happy with a short script that uses a number of tools, e.g. `seqkit` `samtools` `diff` etc.

# Answer
Try md5fa from samtools. It computes two MD5s with and without ordering, respectively. Alternatively, there is seqkit sum, which is perhaps more powerful.
> 3 votes

# Answer
There are a lot of ambiguities in the question, but I'll provide one approach by means of a fasta/fastq sorting program (see here), using process substitution to avoid intermediary files. This will only work for small-ish fasta files, which can fit into computer memory:

```
diff -q <(fastx-sort.pl a.fasta) <(fastx-sort.pl b.fasta)

```

Note that this will compare the entire file. If there are any sequences in one file but not the other, then the result will be that the files are different.

The `-q` option for diff will just output whether or not the two files differ; you didn't specify in your question whether you wanted any more detail. This option can be removed to show differences, or replaced with `-u` to add context. Context and/or line changes are not necessarily going to be helpful for multi-line fasta sequences, because the compared files may have a different order from the input, and the context lines may not include the name of the sequence. This could be fixed by changing the files to single-line fasta format (e.g. by piping through `fasta_formatter` from the FASTX-Toolkit):

```
diff -u <(fastx-sort.pl a.fasta | fasta_formatter) <(fastx-sort.pl b.fasta | fasta_formatter)

```

Note that this approach will output the entire sequence for any different sequences, together with adjacent context (including the sequence name). There are ways round that as well... but it would be far better if you could make your question more specific so that we don't need to explore the universe when answering.

If you want to ignore case, convert everything to upper (or lower) case:

```
diff -q <(fastx-sort.pl a.fasta | perl -pe 'tr/a-z/A-Z/') <(fastx-sort.pl b.fasta | perl -pe 'tr/a-z/A-Z/')

```
> 2 votes

---
Tags: fasta, samtools, seqkit
---


# Question
Title: deseq2 full and reduced model interpretation
I would like to know how to interpret the output of the formula although i been thorough loads of literature but I'm not confident yet.

So this is my full model I'm running

```
dds <- DESeqDataSetFromMatrix(countData=rpkm_ordered, colData=coldata, design= ~ FAB + Sex + Age + TMB + WBC + BM_percentage)

```

So if I have to interpret the output of this how it would be ?

What I have read is this

Here in the above formula we are trying or controlling the `FAB` variable while we look for variation that is coming from other variables?

Now my confusion is when I made this `dds` object the `FAB` was the variable which I fixed while I look for variation coming from other variable.

Next thing is I would like to run a reduced model.

I tried this

```
test.Age <- as_tibble(results(DESeq(dds, test="LRT", reduced= ~Sex + Age )))

```

Now my question is since I'm using the same `dds` object here where the `FAB` is the variable i fixed and with new formula which is

`reduced= ~Sex + Age` this so how would i interpret this model?

Any suggestion or help would be really appreciated

# Answer
By fitting the full model `~FAB + Sex + Age + TMB + WBC + BM_percentage` you can estimate the effect of FAB (or any of the other variables) on RNA expression after adjusting for differences caused by Sex, Age, TMB, WBC, and BM percent. See this post for exactly what it means to adjust for other variables

When you fit a reduced model (in your case `~Sex + Age`) and then conduct a likelihood ratio test (LRT), you are comparing the fit of the full model (`~FAB + Sex + Age + TMB + WBC + BM_percentage`) to the fit of the reduced model (`~Sex + Age`). If the fit is substantially better, the likelihood of the full model will be larger than the reduced model and the p-value will be significant. This means that adding the extra variables that are in the full model (FAB, TMB, WBC, BM\_percentage) improved the fit of the model to the data because together these variables have some type of relationship with RNA expression even after adjusting for the variables in the reduced model.

The challenge with using a LRT is that it doesn't tell you which of the variables added in the full model are responsible for the improved fit - it just indicates that the full model is better. This type of test is useful when you are trying to make broad statements about several variables having a relationship with RNA expression. From the DESeq2 vignette:

> The LRT is therefore useful for testing multiple terms at once, for example testing 3 or more levels of a factor at once, or all interactions between two variables. The LRT for count data is conceptually similar to an analysis of variance (ANOVA) calculation in linear regression, except that in the case of the Negative Binomial GLM, we use an analysis of deviance (ANODEV), where the deviance captures the difference in likelihood between a full and a reduced model.

Definitely read the DESeq2 vignette for more details. In your case if the LR test is significant, it means that FAB, TMB, WBC, and BM percent combined are associated with RNA expression.

Most people are interested in specific comparisons when performing differential expression analyses with DESeq2. These comparisons are usually like do males have higher gene expression than females (comparing sex differences), do patients with higher TMB have greater expression than patients with lower TMB (comparing effect of TMB on RNA expression). If you want to do a specific comparison, you will likely want to use the Wald test, which calculates p-values for individual coefficients (male vs female coefficient, TMB coefficient etc). The Wald test determines whether a coefficient is equal to 0 or some other null hypothesis (see `lfcThreshold` parameter).

The Wald statistic and accompanying p-values are reported by default when using `results()` and not specifying `test="LRT"`

There are some additional differences between Wald tests vs LRT, including accuracy differences when operating on small sample sizes and computational efficiency issues when operating with large sample sizes. Discussion on Wald vs LRT

**TL;DR** Use LRT to test if several variables or factor levels are associated with RNA expression. Use the Wald test to determine if specific comparisons are significant.
> 2 votes

---
Tags: r, deseq2
---


# Question
Title: snakemake tutorial stuck in basics-step-5 (variant calling)
I am learning to use snakemake. I am following the tutorial (https://snakemake.readthedocs.io/en/stable/tutorial/tutorial.html). In the section called "Basics", steps 1 to 4 worked well, processing one sample at a time, from mapping reads until indexing alignments.

```
snakemake --cores 2 sorted_reads/{A,B,C}.bam.bai

```

However, when I add the rule for joint variant calling (step 5):

```
snakemake --cores 2

```

I get the following error:

```
"WorkflowError:
Target rules may not contain wildcards. Please specify concrete files or a rule without wildcards at the command line, or have a rule without wildcards at the very top of your workflow (e.g. the typical "rule all" which just collects all results you want to generate in the end)."

```

It looks like I need to do this "rule all" step, any idea how to set this up?

This is how the snakemake file looks like:

```
SAMPLES = ["A", "B", "C"]

rule bwa_map:
    input:
        "data/genome.fa",
        "data/samples/{sample}.fastq"
    output:
        "mapped_reads/{sample}.bam"
    shell:
        "bwa mem {input} | samtools view -Sb - > {output}"

rule samtools_sort:
    input:
        "mapped_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam"
    shell:
        "samtools sort -T sorted_reads/{wildcards.sample} "
        "-O bam {input} > {output}"

rule samtools_index:
    input:
        "sorted_reads/{sample}.bam"
    output:
        "sorted_reads/{sample}.bam.bai"
    shell:
        "samtools index {input}"

rule bcftools_call:
    input:
        fa="data/genome.fa",
        bam=expand("sorted_reads/{sample}.bam", sample=SAMPLES),
        bai=expand("sorted_reads/{sample}.bam.bai", sample=SAMPLES)
    output:
        "calls/all.vcf"
    shell:
        "bcftools mpileup -f {input.fa} {input.bam} | "
        "bcftools call -mv - > {output}"

```

# Answer
> 2 votes
Read the snakemake documentation. `all` is a pseudo-rule that accepts a set of one or more final output files as input and does nothing with them.

---
Tags: variant-calling, snakemake
---


# Question
Title: How to interpret a negative beta within a binary GWAS?
I have run a GWAS on a binary trait - "Are you following a carnivore diet?" - so the cases are individuals who follow this diet, and the controls are individuals who do not follow this diet. I have a genome-wide significant SNP, with a negative beta in reference to Allele 2. Does this imply that individuals with Allele 2 are less likely to follow this diet?

# Answer
Negative beta means that the trait under test is negatively correlated with the allele. What this actually means in relation to your specific question will depend on how the trait is coded and how the allele is coded.

It might, for example, be the case that "yes" is coded as 0 and "no" as 1, in which case a "no" answer would be associated with a high allele count for allele 2. You will need to check the algorithm to find that out, or (better) do your own manual / independent checking of alleles to discover the direction of association.

Take care when interpreting the results of these association tests, especially when deviating from statements beyond what the results actually say. You should also pay attention to is meaningless, whereas a beta value of -3.43 means there is a strong correlation between the trait and allele.
> 1 votes

---
Tags: gwas, genetics
---


# Question
Title: Where to find database for breast cancer phenotype?
Where can I find databases for breast cancer mutations (SNPs & indels) phenotype, with chromosome, position, ref, alt, gene, clinical significance.

So far, I have only found the Clinvar database.

I need to extract data for breast cancer, in order to train a machine learning model that will be able to predict the pathogenicity of variants (benign, pathogenic etc).

So I did the extraction of data from clinvar vcf file (chromosome, ref, alt, clinical significance), but the clinvar data for breast cancer is very small.

# Answer
> 1 votes
Do you mean this (a product of the BRCA Challenge project):

> https://brcaexchange.org/

---
Tags: database, cancer
---


# Question
Title: How can I get data of single cell RNA sequence with raw count?
I am using dataset GSE85241, but I can't find the read counts of the dataset. It only provides with `RPKM` values.

How can I find the read counts of a dataset?

# Answer
> 1 votes
If the authors did not provide the raw counts, and I don't think it is compulsory, unfortunately you will need to download the raw data, align and count.

---
Tags: rna-seq, single-cell, phylogenetics
---


# Question
Title: Using QCTOOL v2 to process UK Biobank .bgen files - why so slow?
I’m currently using QCTOOL v2 to process imputed .bgen files from UK Biobank, however they seem to be processing very slowly. Is this normal?

My command is pretty basic; I’m filtering out a list of SNPs and samples:

```
/path_to/qctool \
-g /path_to/ukbXXXX_c22_b0_v3.bgen \
-s /path_to/ukbXXXX_c22_b0_v3.sample \
-og /output_path/ukb_c22_filt.bgen \
-os /output_path/ukb_c22_filt.sample \
-excl-rsids /path_to/snps_rem_c22.txt \
-excl-samples /path_to/samples_to_rem.txt

```

It is currently processing SNPs at a rate of **1.2/s** (e.g. 71169/?,57205.8s,1.2/s).

The computational facility I'm using should not limit the speed of an operation.

* Have I made any mistakes in my qctool query?
* Or any other ideas on why it's running so slowly?

Alternative suggestions on how to process these files would be appreciated for example, although I would prefer to use QCTOOL perhaps I have to use PLINK.

---

From the website, https://www.well.ox.ac.uk/~gav/qctool\_v2/

"QCTOOL is a command-line utility program for manipulation and quality control of gwas datasets and other genome-wide data"

# Answer
I contacted the QCTOOL team directly and received this response:

> The design of BGEN means that when you subset samples the data has to be recompressed - this is essentially what makes this slow. (By contrast you can subset SNPs very quickly without recompression using bgenix https://code.enkre.net/bgen.) It is therefore definitely worth considering not subsetting samples but using a sample inclusion/exclusion list instead, if your analysis software supports that.
>
> If you have to subset and want to use QCTOOL - some things to try are: I typically use the options -bgen-compression zstd -bgen-bits 8 now, c.f. https://doi.org/10.1101/308296 this is faster but first check your downstream software supports zstd compression. Use a map/reduce type pipeline (i.e. chunk data for re-encoding) - this can be implemented using bgenix and cat-bgen.

> I assume you’re working with the latest version?
>
> Have you tried first stripping the SNPs, then stripping the samples?
>
> Ever since working with the imputed data from the UKB I was never successful in chopping pieces of bgen files using qctool in a timely manner. If my memory serves me well, qctool manages to strip SNPs fairly quickly, but is super slow to remove samples. Try PLINK 1.9 or 2.0.
>
> As for “The computational facility I'm using should not limit the speed of an operation.”, depending on how disk data moves in/out of slave nodes and how busy the cluster is, I did saw processes to become I/O starved on large clusters.
>
> Also, consider that most tools that perform association testing can utilize SNPs and samples lists. BOLT, SAIGE, regenie, SNPtest… so there might not be a need to pre-filter. Put “NA” as phenotypes you’d like to be “removed” from the testing.

In the end I did not pre-filter SNPs or samples - I set samples to NA within the phenotype file, and used a SNP inclusion list in the second stage of SAIGE with the flag `idstoIncludeFile`.
> 1 votes

---
Tags: plink, qc, ld
---


title: Filter bam with sambamba

text:
# Question
Title: Filter bam with sambamba
I have a very large bam file and I want to filter it to keep only a handfull of positions as defined in a bed file. Can I do this with sambamba or do i need another tool?

Thanks

# Answer
I think it is possible:

```
sambamba view --regions regions.bed --with-header --format=bam reads.bam  > filtered.bam

```
> 1 votes

# Answer
You need `bedtools`, specifically bedtools intersect command, not `sambamba`:

```
bedtools intersect -abam A.bam -b B.bed > A_in_N.bam

```
> 0 votes

---
Tags: sam, sambamba
---


# Question
Title: how to calculate the expression values per gene for all cells
I'm currently working with a seurat object and I'd like to calculate the expression values per gene for all cells within a particular cluster.

I've gone ahead and subsetted the cluster of interest.

```
tumor <- subset(x = data, idents = c("Tumor"), invert = FALSE)

```

How would I calculate the expression values per gene per cell?

```
FindAllMarkers

```

Finds markers (differentially expressed genes) for each of the identity classes in a dataset per ident.

# Answer
> 2 votes
No Seurat function is needed. In the seurat object "cells", you can access the raw counts or normalized counts using `cells@assays$RNA@counts` (raw counts) or `cells@assays$SCT@data` (SCT normalized data)

---
Tags: r, scrnaseq, seurat, gene-expression
---


# Question
Title: Unable to execute DSSP algorithm from BioPython
I have the following settings in PyCharm:

BioPython is properly installed:

**dssp\_test.py**

```
from Bio.PDB import PDBParser
from Bio.PDB.DSSP import DSSP
p = PDBParser()
structure = p.get_structure("16PK", "16pk.pdb")
model = structure[0]
dssp = DSSP(model, "16pk.pdb")
a_key = list(dssp.keys())[2]
print(dssp[a_key])

```

**Output**

```
C:\ProgramData\Miniconda3\python.exe C:/Users/pc/source/repos/DSSP_algoritm_test/dssp_test.py
Traceback (most recent call last):
  File "C:\Users\pc\source\repos\DSSP_algoritm_test\dssp_test.py", line 6, in <module>
    dssp = DSSP(model, "16pk.pdb")
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 437, in __init__
    dssp_dict, dssp_keys = dssp_dict_from_pdb_file(in_file, dssp)
  File "C:\ProgramData\Miniconda3\lib\site-packages\Bio\PDB\DSSP.py", line 231, in dssp_dict_from_pdb_file
    p = subprocess.Popen(
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 947, in __init__
    self._execute_child(args, executable, preexec_fn, close_fds,
  File "C:\ProgramData\Miniconda3\lib\subprocess.py", line 1416, in _execute_child
    hp, ht, pid, tid = _winapi.CreateProcess(executable, args,
FileNotFoundError: [WinError 2] The system cannot find the file specified

Process finished with exit code 1

```

Why am I unable to execute DSSP algorithm from BioPython?

# Answer
> 1 votes
This looks like a python problem and not a “biopython”problem.

Check your file structures. Try giving a full path for 16pk.pdb

These might help. https://stackoverflow.com/questions/17658856/why-am-i-getting-a-filenotfounderror

https://stackoverflow.com/questions/18757127/windowserror-error-2-the-system-cannot-find-the-file-specified-is-not-resol

# Answer
> 1 votes
Ok so this might not provide a solution, but here's a couple things someone might find useful;

* First, Biopython does not automatically contain DSSP, which you have to download as mentioned in the comments below (from here swift.cmbi.umcn.nl/gv/dssp)
* Secondly, DSSP within Biopython, for whatever reason, accepts the path to the file, not the file itself. I see you have that correctly (`dssp = DSSP(model, "16pk.pdb")`), but just pointing it ou
* Third, DSSP originally had `dssp` executable, but the newer version has `mkdssp` executable, so you might need to do `dssp = DSSP(model, "16pk.pdb", dssp="mkdssp")`, so it knows what to associate it with. This was the problem I had, because without specifying the `dssp`, it is not able to launch DSSP.

---
Tags: python, secondary-structure
---


# Question
Title: How to get coding sequences (CDS) of UniProt reference proteomes?
I downloaded the following reference proteomes:

https://ftp.ebi.ac.uk/pub/databases/reference\_proteomes/QfO/

```
'Drosophila':'UP000000803_7227.fasta',
'Arabidopsis':'UP000006548_3702.fasta',
'C.elegans':'UP000001940_6239.fasta',
'Zebrafish':'UP000000437_7955.fasta',
'human':'UP000005640_9606.fasta',
'mouse':'UP000000589_10090.fasta',
'yeast':'UP000002311_559292.fasta'

```

Now, I would like to have the corresponding coding sequences (CDS). My temporary solution is to download CDS files from ENSEMBL FTP server, translate those and match with the ref proteome. However, this only covers 50-90% of protein sequences.

Thanks!

# Answer
> 1 votes
I went through diverse README files and found the answer:

```
Coding DNA Sequence FASTA files (*_DNA.fasta)
=============================================

These files contain the coding DNA sequences (CDS) for the protein sequences
where it was possible.
The format is as in the following example (UP000005640_9606_DNA.fasta):

>sp|A0A183|ENSP00000411070
ATGTCACAGCAGAAGCAGCAATCTTGGAAGCCTCCAAATGTTCCCAAATGCTCCCCTCCC
CAAAGATCAAACCCCTGCCTAGCTCCCTACTCGACTCCTTGTGGTGCTCCCCATTCAGAA
GGTTGTCATTCCAGTTCCCAAAGGCCTGAGGTTCAGAAGCCTAGGAGGGCTCGTCAAAAG
CTGCGCTGCCTAAGTAGGGGCACAACCTACCACTGCAAAGAGGAAGAGTGTGAAGGCGAC
TGA

```

So, EBI itself provides the corresponding coding sequences, which they name **\_DNA.fasta**. I misinterpreted the file-name as genomic (unspliced) DNA.

---
Tags: data-download, sequence
---


# Question
Title: Problem with setting fixed topology for MrBayes
I need to use partitioned DNA input (concatenated multiple genes) and fixed topology (species tree from ASTRAL) to get branch lengths in MrBayes. I have input NEXUS like:

```
#NEXUS
BEGIN TAXA;
  DIMENSIONS NTAX = 318;
  TAXLABELS
    Achillea_nobilis_SAMN11585351
    Anthemis_arvensis_SAMN11585352
    ...
    Haplopappus_sp_SAMN11585375
  ;
END
BEGIN TREES;
  TRANSLATE
    1       Achillea_nobilis_SAMN11585351,
    2       Anthemis_arvensis_SAMN11585352,
    ...
    318     Haplopappus_sp_SAMN11585375
  ;
  TREE treeshrink_exons_good_nobrl = [&R] (1,((2,3...);
END;
BEGIN DATA;
dimensions ntax=318 nchar=169024;
format datatype=dna missing=? gap=-;
MATRIX
Achillea_nobilis_SAMN11585351               TTATGCAA...
Anthemis_arvensis_SAMN11585352              TTATGCAA...
...
Haplopappus_sp_SAMN11585375                 TTATGCAA...
;
END;
BEGIN MRBAYES;
  charset At1g09700 = 1-213;
  charset At1g10740 = 214-456;
  ...
  charset At5g66280 = 168530-168757;
  charset At5g66680 = 168758-169024;
  outgroup Achillea_nobilis_SAMN11585351;
  partition by_gene=467:At1g09700,At1g10740,...,At5g66280,At5g66680;
  set partition=by_gene;
  prset topologypr=fixed(treeshrink_exons_good_nobrl);
  prset applyto=(all) ratepr=variable;
  lset applyto=(all) nst=mixed rates=gamma ngammacat=10;
  unlink tratio=(all) revmat=(all) omega=(all) statefreq=(all) shape=(all) pinvar=(all) correlation=(all) ratemultiplier=(all);
  mcmcp ngen=10000000 nruns=4 nchains=4 samplefreq=5000;
  ss;
  sump relburnin=yes burninfrac=0.25 nruns=4;
  sumt relburnin=yes burninfrac=0.25 nruns=4;
END;

```

I'm not 100% sure about the `MRBAYES` block, but generally I think it's correct. When I launch it (MrBayes 3.2.7a) I get an error which I'm unable to resolve:

```
$ mb sp_treeshrink_exons_good.mb.nex
...
   Executing file "sp_treeshrink_exons_good.mb.nex"
   UNIX line termination
   Longest line length = 169069
   Parsing file
   Expecting NEXUS formatted file
   Reading taxa block
      Allocated taxon set
      Defining new set of 318 taxa
   Exiting taxa block
   Reading trees block
      Successfully read tree 'treeshrink_exons_good_nobrl' # NOTE
   Exiting trees block
   Deleting user trees
   Deleting previously defined taxa
   Reading data block
      Allocated taxon set
      Allocated matrix
      Defining new matrix with 318 taxa and 169024 characters
      Data is Dna
      Missing data coded as ?
      Gaps coded as -
      Taxon   1 -> Achillea_nobilis_SAMN11585351
      Taxon   2 -> Anthemis_arvensis_SAMN11585352
...
      Taxon 318 -> Haplopappus_sp_SAMN11585375
      Successfully read matrix
      Setting default partition (does not divide up characters)
      Setting model defaults
      Seed (for generating default start values) = 1654610114
      Setting output file names to "sp_treeshrink_exons_good.mb.nex.run<i>.<p|t>"
   Exiting data block
   Reading mrbayes block
      Defining charset called 'At1g09700'
      Defining charset called 'At1g10740'
...
      Defining charset called 'At5g66280'
      Defining charset called 'At5g66680'
      Defining partition called 'by_gene'
      Setting by_gene as the partition, dividing characters into 467 parts.
      Setting model defaults
      Seed (for generating default start values) = 1901076251
      Setting Topologypr to Fixed for partition 1
      Setting Topologypr to Fixed for partition 2
...
      Setting Topologypr to Fixed for partition 466
      Setting Topologypr to Fixed for partition 467
      Could not find tree 'treeshrink_exons_good_nobrl' # NOTE
      Could not set fixed topology from user tree 'treeshrink_exons_good_nobrl' # NOTE
      Error when setting parameter "Topologypr" (2)
      The error occurred when reading char. 26-52 on line 1445
         in the file 'sp_treeshrink_exons_good.mb.nex'
...

```

So it correctly reads tree `treeshrink_exons_good_nobrl` and after while it complains it couldn't find tree `treeshrink_exons_good_nobrl`...? I can't find where this weird error comes from.

# Answer
The answer is pretty simple. What you are requesting is fixing the tree prior at a tree you are specifying. The getter/setter is complaining it cannot fit the data to the specified tree because of 26 bp at 26-52 on line 1445 of the nexus file.

The first thing needed is to post output to:

```
sed -n '1445p' ./sp_treeshrink_exons_good.mb.nex

```

What I suspect you will find are non-nucleotide characters associated with "missing data" for example `xxxxx`, or `0000` what the nexus file wants is `?????` for missing data. It is possible part of the degenerate nucleotide code is not being observed (long while since I used MB), e.g. `RYB` I think its unlikely.

The second approach, if the output of the `sed` command is pretty normal is to remove the prior setting the tree, i.e. the tree you are supplying is violating e.g. Jukes-Cantor correction.

Finally and a related issue is there's an alignment error, i.e. the alignment at this region is out of frame for these residues and Jukes-Cantor is violated so it refuses to proceed. This would be really unusual however, because I've used MB where there are clear pockets of saturation and it hasn't complained.

Personally, I'm impressed by the level of foresight in MB under written into this getter/setter (its checking your input). This algorithm is rarely used but is clearly well written.

---

Okay there is no line 1445. What I recommend is placing the alignment in either `Sequtron` or `Clustal X` and scrolling down at the position 26-52 and if not successful 'optical analysis' of the alignment.

Personally, I would use Biopython AlignIO for this task and `Counter` from the `Collections` library to interrogate the bug. However, thats coding stuff and optical analysis does the same thing. Essentially, you loop through each sequence and `Counter` will make a frequency count of everything in each sequence. You then say if you are not `AGCT-` tell me what taxa this is. If you code in Python thats a good follow-up question (separate question).
> 0 votes

---
Tags: phylogenetics, phylogeny, mrbayes
---


# Question
Title: Should I expect / use a prokaryote-specific substitution model when building a phylogenetic tree?
I have a couple of prokaryotic genes that I'm trying to build phylogenetic trees from.

I inputted them into IQTree, which, as part of the workflow, runs through multiple substitution models until it finds the one that seems to fit the data best.

My issue: quite a few of the genes ended up with eukaryotic\_specific substitution models (such as Q.plant or Q.insect). Is this likely to bias my data? Should I re-run the analysis and remove the substitution models that are meant for different groups?

Thank you for your time!

# Answer
> 3 votes
IQTree is a fashionable method for speeding up bootstrapping and used on ultra large datasets. What IQTree might be doing is rigidly fixing a model and parameterisation is not estimate via the ML. This might be how it achieves its performance.

In context, IQTree can use an establised mixed model, which is one approach to protein phylogeny. As a personal note, I don't favour it and come back to the reason later. It MIGHT also generate its own mixed model. There are not 700 established (published) amino acid matrices plus established mixed models, combined there are no more than 50.

There are around 30 protein mutation matrices, Dayhoff models (I honestly don't who Dayhoff is), one for insects which might be for mtDNA ... mtArt, one for mtDNA, mtREV I think and so on. Please check carefully the number of models it is evaluating.

Its easy to see how 4 or 5 combinations (or possibly permutations) of 30 become 700 models if it is generating its own mixed model.

Alternatively, what it does (mixed model or not) BIC and AIC do is given a tree estimate (prior) what's model is the best probabilistic fit. In your case Q.insect is cropping up, whether this is mixed model or not ... dunno.

What I don't know is how it calculates the mixed model, is it creating a new singular matrix by averaging (which is probably closer to the algorithm), or is it sequentially passing the data through first one model, then another etc ... Personally, I avoid mixed models for precisely this reason, I dunno what the actual matrix is, but accept it can provide strongest AIC/BIC.

It suggests either:

* BIC/AIC results in an outcome that isn't biologically realistic .... **OR**
* there has been some weird parallel/convergent evolution - this happens in parasites for certain

If you are using a single matrix not mixed model and Q.insect is cropping up - **I would be concerned**. If it is part of a mixed model, I wouldn't worry.

Is it checking 30 models (plus published mixed models) or 700 matrices??

If the final point was correct (biologically real result), you check it via blast - simply exclude bacteria and see if you get a plant of insect sequence as the top hit.

Ultimately, the uncertainty is the reason I avoid mixed models.

---

To answer the question, I'm not aware of a bacterial amino acid matrix. HOWEVER, you can select a mixed model. It might just be screening around 50 models via AIC/BIC and within that Q.insect is the most probable.

**Summary** The alternative is to look for support using a general amino acid matrix (for example JTT, LG or JTTDCMut).

Personally, I would be very cautious about proceeding with a matrix generated by insect genomes outside a mixed model for bacteria and in the absence of any biological support. FOR EXAMPLE, if these were insect symbionts - that makes sense.

Which ever approach is used a Bayesian approach is often required.

---
Tags: phylogenetics, phylogeny, bacteria, substitution-model, matrix
---