# Question
Title: Running cellranger on scRNASeq data with feature barcoding (10x + antibody capture)
I can't seem to find a clear answer to this question, so here it goes:

I have sequenced scRNASeq + scVDJSeq (TCR) data, which has been sequenced using feature barcoding from 10x genomics, via antibody capture.

There is a very useful `cellranger multi` command, which allows us to process both this data types in 'one go'.

In order to run it one needs to generate a `csv` file a priori, such that `cellranger` is able to identify which files are pertinent to which file type - see this link:https://support.10xgenomics.com/single-cell-vdj/software/pipelines/latest/using/multi

In the link above, with the example given under "Example multi config CSVs", it seems as though there are separate `fastq` files for the feature barcode, versus for the actual gene expression samples. Whereas, in my case, I only have the standard `R1` and `R2` fastq files, which is creating some confusion as how to generate this file, and run cellranger.

Now my question is, how does one run `cellranger multi` with feature barcoding should there only be the standard 2 `fastq` files? My understanding is that the feature barcode is usually in the `R2` read? Is this correct? Would I need to Demultiplex? I am completely lost as to how to proceed. If you know, could you please provide a short example of what your `csv` files would look like, should you not have separate `fastq` files for the feature barcode files?

Thanks,

# Answer
Check your FASTQ headers - if they've been demultiplexed, the index barcodes can be found in the header. The barcodes are usually separated by a `+` in the comment part (second segment when the header is split by `" "` blank space) of the header. If you see a match to `[ATGC]{10}[+][ATGC]{10}` in the header, the first `[ATGC]{10}` is the I1 barcode and the second is the I2 barcode. You can use `FFFFFFFFFF` placeholder qual scores to generate "fake" I1 and I2 files and then use that with cellranger.
> 2 votes

---
Tags: scrnaseq, fastq, 10x-genomics, cellranger, demultiplexing
---


# Question
Title: Help with analyzing SNP data for three patient groups
I am currently working on a project where I have genetic data for patients classified into three different groups. Specifically, I have information on 26 SNPs for each patient. I have recorded each allele using letters. I am trying to decide on the best approach to analyze this data and would appreciate any guidance or suggestions. Firstly, I need help formatting the SNP information for analysis. I have recorded each allele using a comma-separated format (e.g., "A,G"). However, other encoding methods may be more appropriate for analyzing SNP data. Could someone please suggest a suitable encoding method for my data? Additionally, I am interested in performing machine learning and cluster analysis on this dataset to investigate the relationship between genotype and patient group. Please recommend any specific algorithms or packages that may be useful for this type of analysis in either Python or R. Thank you in advance for your help and guidance.

# Answer
**Overview** In the first instance unsupervised learning is needed - namely PCA - and strongly recommended to identify a training target. The point of this analysis is the resulting clusters of different 'patients' (?).

**How to do it?** The first part of the analysis is to perform one-hot-key coding.

```
import pandas as pd
df = pd.read_csv({'A':['A', 'G', 'C'], 'B':['T', 'A', 'G']})
one_hot_encode = pd.get_dummies(df, columns = ['A', 'B'])

```

This converts A, G, C, T to numerical values and then gets passed to unsupervised learning.

**Supervised learning** Moving into supervised ML is not trivial, even with a suitable 'target', unless a simple algorithm such as naive Bayes is used. First I'd do the PCA analysis. ...

**Edit** the training target could simply be the three patient groups.
> 0 votes

---
Tags: snp, correlation
---


# Question
Title: Extract specific nucleotide base position for query and subject from text file
I have an input file looks like this:

```
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

```

The input file is the result file after running the command line blast. It's a txt file generated after blast alignment and in txt file header format are like this.

I am trying to extract the specific nucleotide base for position 129, 249, 366, 372 and 392 from query and subject aligned sequence. python and perl solution will be fine as well. thank you.

For example: I need my output to have format below like this:

```
Seq_name        base_position       query_base  subject_base    all_subject_base
CD630_08510_19  129                 G           T               TAT-C
CD630_08510_19  249                 G           A
CD630_08510_19  366                 C           T
CD630_08510_19  372                 A           -
CD630_08510_19  392                 T           C

CD630_08510_20  129                 G           G               GATCC
CD630_08510_20  249                 G           A
CD630_08510_20  366                 C           T
CD630_08510_20  372                 A           C
CD630_08510_20  392                 T           C

```

all\_subject\_base: are the bases from position 129, 249, 366, 372 and 392 of subject\_base column. I want to keep this in one column thats why all\_subject\_base.

I tried to extract this with grep but have not been able to. I am new in bioinformatics so any suggestion would be very helpful.

# Answer
> 2 votes
Usually, when you run your BLAST alignment you specify the output format (e.g. as described in "outfmt" in https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options\_common\_to\_all\_blast/). If, for some reason, you can't re-run the BLAST command (e.g. if this output file was given to you, or if you can't access your query or reference sequence/s) you *can* parse this 'readable' output format but it's not straightforward.

Here is a potential solution using awk:

```
cat positions.txt
129
249
366
372
392

cat input.txt
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
> CD630_08510_20
Query  121  TGGCTTGGATTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-CT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-CT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

```

If you have more positions:

```
cat positions.txt
129
249
366
372
373
374
375
376
392

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-GGTACT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  373 G   G
CD630_08510_19  374 G   G
CD630_08510_19  375 T   T
CD630_08510_19  376 A   A
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-GGTACT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  373 G   G
CD630_08510_20  374 G   G
CD630_08510_20  375 T   T
CD630_08510_20  376 A   A
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

```

NB this won't work if you only have a single position in positions.txt, but if that is a potential use-case you can modify the script to suit, e.g. replace `(j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")` with `f[j]`.

---
Tags: blast, linux, shell
---


# Question
Title: Using ticks on the x axis or a scale bar when visualisizing a phylogram?
When visuliazing a phylogram is it preferable to show ticks on the x axis or a scale bar?

Given an example tree with branch lengths in newick format `(A:0.01,(B:0.02,C:0.05):0.02);`, the two options are illustrated in the plot below:

# Answer
> 1 votes
Both will work. The convention is the scale-bar and the rationale is below (which is a bit complicated).

By convention:

* the right-hand tree where there is a scale-bar and is the traditional diagram for phylograms.
* the left-hand tree scale-bar ("x-axis") is used for molecular clock representations, i.e. not for phylograms, but dendrograms.

In terms of clarity the "x-axis" (or extended scale bar) appears clearer. However, firstly, for a seasoned phylogenecist it makes no difference and simplicity is preferred. Secondly, there's coalescent theory.

The left-hand tree ... the thing is in a molecular clock representation the scale is a **reverse axes** starting at the right-hand tips. Whether thats time or genetic diversity its a good representation (well for time its essential) *because* coalescent theory means that there is zero genetic diversity for that taxa at time 0 (the idea is you are extrapolating, or reconstructing, the common ancestor). The representation you presented confuses that convention because 0 starts at the assumed node in the outgroup.

Thus the scale-bar is better.

---
Tags: phylogenetics, phylogeny, visualization
---


title: Manta "--exome: command not found"

text:
# Question
Title: Manta "--exome: command not found"
I am using Manta (https://github.com/Illumina/manta/blob/master/docs/userGuide/README.md) but whenever I try to specify the `--exome` command using the following code I get the error below:

The code in my script:

```
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir}
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

```

The `runDir` is directly below the previous command however, it doesn't seem to write well in this post.

I ran the script via the following command:

```
./script.sh /path/to/reference.fa /path/to/bam_files /path/to/output_dir

```

I get the error:

```
/path/to/my/sample/folder/sample1/runWorkflow.py
/path/to/this/job_scripts/number_by_HPC: line 40: --exome: command not found

```

How do I get it to run with the `--exome` command and not get the following error (it works without the `--exome` command)? I am using exome data so I believe I need to use the `--exome` command

# Answer
You're missing a `\` in the line before the one containing the `--exome`, so shell treats it as a new command, which is why you see this error. Add the  `\` after `${sub_dir}` and your script should not give you this error.

With that correction added in, your code will look like:

```
# Accept path to reference fasta file as a parameter:

reference=$1

# Accept input and output directories as parameters:

input_dir=$2

output_dir=$3

# Hard coded Manta Installation Path:

## If someone installs Manta they must specify where the installation is placed:

MANTA_INSTALL_PATH="/path/to/manta-1.6.0.release_src/install"

# Main code:
for file in ${input_dir}/*.bam
do
  # Create subdirectory named after bam file
  sub_dir=${output_dir}/$(basename ${file} .bam)
  mkdir -p ${sub_dir}

  # Use reference fasta path specified by user
  ${MANTA_INSTALL_PATH}/bin/configManta.py \
    --bam $file \
    --referenceFasta ${reference} \
--runDir ${sub_dir} \
    --exome
done

# Now run each `runWorkflow.py`

find "$output_dir" -name "runWorkflow.py" -exec {} \;

```
> 4 votes

---
Tags: bash, exome
---


# Question
Title: Passing strings for password, email authentication and download link
Hi I am downloading the data from a website which requires authentication of email\_id and password before the download. They then provide the download link which is used to download dataset. This download link is different every time

```
### Download Cosmic Files

##First you need to login with your login details and the access key ID may vary with each login ID

##STEP-1
echo "pr@gmail.com:abc123@" | base64
#The above lines gives a authorization key cHIxNEBzYW5

##STEP-2
# get the download link using the authorization key from step 1
#curl -H "Authorization: Basic cHIxNEBzYW5" https://cancer.sanger.ac.uk/cosmic/file_download/GRCh38/cosmic/v97/CosmicSample.tsv.gz

#The above line gives the download link where AWS access key changes everytime like in this case it was https://cog.sanger.ac.uk/cosmic/GRCh38/cosmic/v97/CosmicSample.tsv.gz?AWSAccessKeyId=KRV7P7QR9DL41J9EWGA2&Expires=1679591113&Signature=WPjwDQsMoaLtRA2utiGLNrOSlOA%3D

##STEP-3
# Download file with the link from step 2
#curl -o CosmicSample.tsv "https://cog.sanger.ac.uk/cosmic/GRCh38/cosmic/v97/CosmicSample.tsv.gz?AWSAccessKeyId=KRV7P7QR9DL41J9EWGA2&Expires=1679591113&Signature=WPjwDQsMoaLtRA2utiGLNrOSlOA%3D"

```
1. I want to automate the script first by passing used\_id and password as an string
2. I want to replace the AWS key in the step 3 with a string that is generated from step 2

Thank you

# Answer
> 3 votes
The cURL command to obtain the download link returns a JSON containing the URL. This can be parsed on the command line using jq and used to construct a command to download the file. For example, with a list of URLs in a file (one per line) called *urls.txt*, you could try the following:

```
while IFS= read -r url ; do
    >&2 echo "Downloading: ${url}"

    curl -sSH @cosmic_header "${url}" |
    jq --arg base "$(basename "${url}")" -r '"curl -o \($base) \(.url | @sh)"' |
    sh
done < urls.txt

```

Note that the above stores the HTTP header in a file called *cosmic\_header*, which can be created using:

```
echo "Authorization: Basic $(echo 'username:password' | base64)" > cosmic_header

```

This avoids exposing the authorization key on the command line when you go to download your list of files, which might be preferable if working in a shared environment.

---
Tags: database, bash, curl
---


# Question
Title: Improving conversion of abricate tsv file to gff3 file
Since such a neat solution (abricate tsv to gff3) was provided by Steve, here are few other steps that I am looking to add so that the script progress to logical maturity to be usable by many others.

I have two files - (1) fasta file with .fna extension, and (2) a data file with .tsv extension generated from abricate.

I need to write a gff3 file (containing sequences) by combining these two files (.fna and .tsv)

*A .gff3 file generally has three or four components as below.* Details here.

1. Declaration of file type and version
2. Declaration of sequence region(s) starting from Name of the sequence, followed by single space, then 1 separated by single space, ending in numerical value of the sequence length of the fasta sequence (which is the .fna file).
3. Sequence features
4. Fasta sequence (optional)

See below a truncated example

##gff-version 3
##sequence-region NZ\_BHEI01000001 1 13890
NZ\_BHEI01000001 Prodigal CDS 29 1588 . + 0 ID=CEMHHHDG\_00001;product=hypothetical protein
NZ\_BHEI01000001 Prodigal CDS 1585 2289 . + 0 ID=CEMHHHDG\_00002;product=Putative beta-lactamase HcpC
....
##FASTA
\>NZ\_BHEI01000001
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........

Notice that the sequence names has to match in every instance (declaration of sequence regions, sequence features, and the fasta header at the end)

The .tsv to .gff3 conversion was elegantly solved by Steve.

Now we need to add the declaration of sequence region(s) in the .gff3 file. The declaration comprises of
"##sequence-region"\[space\]the input file name without the file extension .fna\[space\]sequence start\[space\]sequence length\[\n\]

For the last part, I have used the following awk code on the .fna file to generate the length in desired format.

```
awk '/^>/ {if (seqlen); seqlen=0;next; } { seqlen += length($0)}END{print "1 " seqlen}' example.fna

```

This output needs to be appended in the declaration of sequence region line.

Now we need to append the fasta sequence in the following format

##FASTA
\>Sequence name as in the parts above in the .gff3 file
AAAAAAAGAATGCTTGCCTACTGGAGTAGTGCCTGTTACTTTTAAACCCACTTTTTTGCG
TTATGCGGATGAGTATTTTTTAAGAGTTGAATTTCAAGATGGAAGTGATGAAATCACTCA
TATAGAGGAGTTGGCAAAATACACAGATCAAA..........

**Is there any help in sight, please?**

# Answer
> 1 votes
ABRicate can handle multiple input files and will set the *FILE* column in it's output accordingly. Assuming all of the FASTA files in this column exist on disk, you can try the following solution. It uses pyfaidx, which allows for fast random access to the sequences in the (indexed) FASTA file(s). If the FASTA indexes don't exist, they will be created. It loads only a minimal amount of the input FASTA file(s) into memory. With the following in a file called `tsv2gff3.py`:

```
import argparse
import csv
import contextlib
import re
import sys
import textwrap

from dataclasses import dataclass, astuple
from pathlib import Path
from typing import List, Literal, Tuple

from pyfaidx import Fasta

```

```
def get_argument_parser():
    parser = argparse.ArgumentParser(add_help=False)
    group = parser.add_argument_group('positional arguments')
    group.add_argument(
        "tsv",
        type=str,
        metavar="FILE",
        help="The ABRicate TSV file",
    )
    group = parser.add_argument_group('output options')
    group.add_argument(
        "-o",
        "--output",
        type=str,
        metavar="FILE",
        default='-',
        help="Write the GFF3 output to FILE (default: stdout)",
    )
    group.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )
    group = parser.add_argument_group('additional options')
    group.add_argument(
        "-h",
        "--help",
        action="help",
        help="Show this help message and exit",
    )

    return parser

```

```
@dataclass
class Feature:

    seq_id: str
    source: str
    type: str
    start: int
    end: int
    score: float
    strand: Literal['-', '+']
    phase: Literal['0', '1', '2']
    attrs: List[Tuple[str, str]]

    @property
    def attributes(self):
        return ';'.join([f'{key}={val}' for key, val in self.attrs])

    def __iter__(self):
        return iter(astuple(self))

    def __str__(self):
        vals = list(map(str, astuple(self)[:-1])) + [self.attributes]
        return '\t'.join(vals) + '\n'

```

```
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:
        if args.output != '-':
            out = stack.enter_context(open(args.output, 'wt' if args.force else 'xt'))
        else:
            out = sys.stdout

        tsv = stack.enter_context(open(args.tsv, 'r'))
        first_line = tsv.readline().rstrip()
        fieldnames = first_line.lstrip('#').split('\t')
        reader = csv.DictReader(tsv, dialect='excel-tab', fieldnames=fieldnames)

        fasta_files = dict()
        seq_ids = dict()

```

```
        out.write('##gff-version 3\n')
        for idx, row in enumerate(reader, 1):

            filename = Path(row['FILE'])
            if filename not in fasta_files:
                fasta_files[str(filename)] = stack.enter_context(Fasta(filename))

            if filename.suffix in ('.fasta', '.fa', '.fna'):
                seq_id = filename.stem
            else:
                seq_id = filename.name

            fasta = fasta_files[str(filename)]
            seq = fasta[row['SEQUENCE']]

            if not seq_id in seq_ids:
                out.write(f'##sequence-region {seq_id} 1 {len(seq)}\n')
                seq_ids[seq_id] = [row['FILE'], row['SEQUENCE']]

            database = row['DATABASE']
            product = re.sub(r'[\[\]\(\)]', '', row['PRODUCT'])
            start, end = int(row['START']), int(row['END'])
            attrs = [
                ('ID', f'{seq_id}_{idx}'),
                ('Product', product),
            ],

            feature = Feature(seq_id, database, "CDS", start, end, '.', '+', '0', *attrs)

            out.write(str(feature))

```

```
        out.write('##FASTA\n')
        for seq_id, (filename, seq_name) in seq_ids.items():
            fasta_rec = fasta_files[filename][seq_name]
            seq = textwrap.fill(str(fasta_rec), width=60)

            out.write(f'>{seq_id}\n{seq}\n')

```

```
if __name__ == '__main__':
    main()

```

And `pyfaidx` installed, the usage is:

```
$ python tsv2gff3.py -h
usage: tsv2gff3.py [-o FILE] [-f] [-h] FILE

positional arguments:
  FILE                  The ABRicate TSV file

output options:
  -o FILE, --output FILE
                        Write the GFF3 output to FILE (default: stdout)
  -f, --force           Overwrite the output file if it exists

additional options:
  -h, --help            Show this help message and exit

```

---
Tags: gff3, tsv
---


# Question
Title: Analysing the Effect of a Drug on the Morphological Changes of a Cell Type
We came across a project in our lab that no one exactly knows how to approach. Since, I know a little bit of Python programming, this project was assigned to me.

There is a data from a randomised controlled clinical trial with 60 participants. Half (30) are control and the other half are actual patients (case). From each group, biopsy samples were taken before (pre) and after (post) treatment with drug A that targets cell type B in the tissue.

In each biopsy sample, the B cells were detected and classified either as “normal” or “malignant”.

Associated metadata for each biopsy are displayed here (I have only included 12 records for this file):

```
name    patient_id  arm treatment
111     0       control     pre
112     1       control     pre
113     2       control     pre
121     0       control     post
122     1       control     post
123     2       control     post
211     75      case        pre
212     76      case        pre
213     77      case        pre
221     75      case        post
222     76      case        post
223     77      case        post

```

● name: spreadsheet file name

● patient\_id: patient identity number

● arm: trial arm (‘case’ or ‘control)

● treatment: treatment condition (‘pre’ or ‘post’ treatment)

Other files (I have only included 20 records for each file: 10 normal and 10 malignant) contain cell detection results from a single biopsy, and each row in a spreadsheet represents an individual detected cell:

● x: x coordinate

● y: y coordinate

● label: classified label (‘normal’ or ‘malignant’)

Just showing an example, file 111 looks like this:

```
x   y   label
730 724 normal
1962 450 normal
1511 817 normal
1244 455 normal
2529 397 normal
1878 262 normal
2248 369 normal
2007 273 normal
1531 878 normal
1729 834 normal
931  1270 malignant
1282 314 malignant
1630 839 malignant
1543 460 malignant
2493 237 malignant
1311 744 malignant
1999 366 malignant
737 1361 malignant
2252 448 malignant
2620 398 malignant

```

The rest can be found here, but probably they will not be necessary: https://www.mediafire.com/file/ka7r59kf0swnbnd/OtherFiles.rar/file

I am trying to answer 3 questions here (if you can think of other questions, please let me know):

1. identifying post-treatment morphological changes due to the effect of drug A.
2. proposing measures to quantify the changes.
3. Using appropriate statistical analysis tools to provide insight whether the changes are due to chance or not.

# Answer
> 1 votes
Yes it's easy you load everything into a dataframe, so each variable is one column and perform supervised machine learning. That will solve all questions.

**Hypothesis 1** The target is 'pre' 'post' **if there is only one drug in the trial**. If there is more than one drug in the trial the drugs would become the target, e.g. if you wanted to know the efficacy between 'drug A' and 'drug B'. What you want to know is the effective of the drug. This would help hypothesis 2.

**Hypothesis 2** Estimated cell size is the training target. This is more complicated because I'm not sure what x and y mean: but if x and y are cell measurements then an estimate of cell surface area should be made. If x and y can become a single variate of cell surface area than that can form categorical data ... super size, big, medium large, medium, medium-regular, regular etc ... This could then become the training target and probably what you are looking for.

The ML feature selection will assess cell size (if 'pre' and 'post' is the training target; or), or 'pre'-'post' (if 'cell size estimation' is the training target). In either case the expected result is it should provide the strongest weight in the analysis.

**Finally** Finally and importantly, the more variables you can input the better, e.g. the gender of the patient, age of the patient as categorical data (bins), smoking, drinking ... anything you can think of, the more the better. ML should ignore anything that has not impact on 'pre' / 'post' or 'estimated cell size'.

The ML algorithm is more complicated but start with naive Bayes and see how it goes. Beyond this it does start to get complicated however.

---

To address the comments,

* Load everything into a dataframe

A dataframe is a basic unit of operation in `Python` (`pandas`) or `R`. Python also uses `numpy` arrays. Everything, means everything - all the data you have, metadata all files. For the calculation the columns involving patient id and patient number, this might be better excluded. I'm not sure - it's complicated because you'll have duplication of records - I would do the calculations with and without this data. You can delete patient id (it's the same as patient number - from a calculation point of view). The x, y is a special case.

Thus all the various tables you describe get loaded into one dataframe, this includes all the patient data such as gender, weight, age, social factors ... it's just like one giant excel table, with rows and columns. The columns are labelled according to the variable e.g. Drug therapy = 'pre' and 'post' data. The dataframe will feed directly to the ML algorithm - a good reason for using it.

***x*** and ***y***

Personally I don't think splitting this data into two columns is a good idea. It will confuse the algorithms - its needs to be one column of data - a single variate. To achieve this will be doable. I would use **π** to calculate the cell surface area. How that is done for an elongated circle - I don't know - but someone will have done it somewhere, you just need to look it up. Thus 2**π***r*2 if its perfectly round where *r* = diameter/2 ... elongation must have a derivation its just finding it.

This will absolutely make a difference in your calculations - what I would do is make it to a categorical variable by plotting it and constructing a histogram - otherwise you are into transformations and thats complicated. Categorical variables will also permit logistic regression ML and that can be very useful, because it places each variable into positive or negative weights.

---
Tags: python, drugs, meta-analysis
---


title: Dealing With Manta Limitations

text:
# Question
Title: Dealing With Manta Limitations
I am trying to figure out a way to remove the limiting factors that Manta says it cannot handle. To quote from the Manta page:

> The following limitations exist on the input BAM or CRAM files provided to Manta:
>
> * Alignments cannot have an unknown read sequence (SEQ="\*")
> * Alignments cannot contain the "=" character in the SEQ field.
> * Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.
> * Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)

Would using a `find` statement in this way and then making `.bai` files be adequate to remove the limiting factors that Manta has?:

```
find ${input_dir} -name "*.bam" -type f -exec sh -c 'output_file="${0%.bam}_filtered.bam"; temp_file=$(mktemp); samtools view -h -q 70 "$0" | sed "/^\(@\|[^[:blank:]]\+[[:blank:]]\+4\)/!s/\t[=*]/\t/g" | samtools view -Sb - > "$temp_file"; mv "$temp_file" "$output_file"; samtools index "$output_file"; rm "$0"' {} \;
# Make the .bai files:

for file in ${input_dir}/*.bam; do
    samtools index "${file}" "${file%.*}.bai"
done

```

However, I get the error (repeatedly for all my bam files):

> \[main\_samview\] fail to read the header from "-". samtools index: "/path/to/MyBatch1And2/SampleNumber\_filtered\_filtered\_filtered.bam" is in a format that cannot be usefully indexed

I am not sure why I got this error once and then it did not occur again. Why? Or would this not be enough remove the limiting factors (as there are more to be considered or the code is not written correctly) of Manta?

# Answer
Assuming you have **paired-end** Illumina reads mapped using BWA or BWA-MEM (with default parameters), these limitations are usually not much of a concern:

> * Alignments cannot have an unknown read sequence (SEQ="\*")

Best practice for representing alignments in SAM format, suggests that the SEQ and QUAL of secondary alignments should be set to an asterisk to reduce file size. You can either filter these out after the fact, using something like `samtools view -F 256 <bam>`, or simply avoid writing out any secondary alignments in the first place.

---

> * Alignments cannot contain the "=" character in the SEQ field.

An "=" character can be used to denote that the base is identical to the reference base. I've not come across an aligner or downstream process that produces such output. If you know of such a tool or process, please comment below to let me know.

---

> * Alignments cannot use the sequence match/mismatch ("="/"X") CIGAR notation

Please see: Why most aligners do not output the "X" CIGAR operation?

---

> * RG (read group) tags in the alignment records are ignored -- each file will be treated as representing one sample.

If you have multiple libraries for a single sample, and they differ in some not so small way (perhaps you have a paired-end library and an old mate-pair library, or maybe the distribution of insert sizes/template lengths of each library varies wildly), you may want to consider running each library through separately.

---

> * Alignments with basecall quality values greater than 70 are rejected (these are not supported on the assumption that this indicates an offset error)

Short-read sequencing technologies do not produce such high base quality scores. Only long-read technologies, like PacBio and Oxford Nanopore, produce base quality scores above 70 (up to 93).

---

---

Ultimately, I don't think you really need to deal with any of these with some command pipeline if you choose your inputs wisely.
> 2 votes

---
Tags: sam, bash
---


# Question
Title: RNAseq alignment: best practices for aligning to multiple isoforms?
I have Illumina RNAseq data and would like to maximize my power to find candidate genes that are differentially expressed genes between experimental conditions.

Many of my (de novo assembled and annotated) gene models have multiple isoforms. Because I will only be using "uniquely" aligning read counts (as specified by DESeq2 documentation) I fear that I might be disadvantaging myself by aligning my reads back to the full set of gene annotations.

Would I be wrong to first filter my gff file for only the longest isoform prior to read alignment?

# Answer
The fact that an alignment is unique or not is relative to the genome, not its annotation. In other words, a read will be considered a multimapper (and ignored by some programs) if it can align to several genomic positions, but a read will be considered properly mapped if it aligns to a single genomic position, that may or may not correspond to multiple isoforms.

So if you're only interested in gene expression, there shouldn't be any problem with multimappers.

Note that this answer could change depending on what you are aligning on, and how exactly you performed your *de novo* assembly.
> 2 votes

# Answer
If you want to account for changing abundances of isoforms, use a pseudoaligner like kallisto for counting. And use DEXSeq, not DESeq2, as DESeq2 is intended to be applied to gene counts, not transcript counts.
> 1 votes

---
Tags: rna-seq, sequence-alignment, deseq2, differential-expression, isoform
---


# Question
Title: What is the NCBI's definition of an "atypical genome"?
Using the new NCBI Datasets platform, you can browse the collection of genomes associated with one or more taxa. For example, searching *Pseudomonas aeruginosa* returns 19,878 genomes as of 29 March 2023.

In the search filtering tab, they give the option to "Exclude atypical genomes". Applying this option to the previous search returns 19,737 genomes, implying that 141 genomes were discarded as atypical.

**What criteria are they using to label genomes "atypical"?** I've read the documentation for the NCBI Prokaryotic Genome Annotation Process, and my hunch is that atypical genomes are those that are annotated with atypical genes, though I'm not sure.

# Answer
> 3 votes
Biostars user MirianT\_NCBI provided an answer on my crosspost, which I'm copying here:

> You can find more information about what constitute an atypical genome in this FAQ. Here's the list of reasons for a genome to be marked as atypical:
> Datasets flags genomes with a warning icon and message for the following genome problems:
>
> * chimeric - sequences from two different organisms are joined together.
> * contaminated - sequences from another organism, cloning vectors, linkers, adapters or primers are present in the assembly.
> * genome length too large - total non-gapped sequence length of the assembly is more than 1.5 times that of the average for the genome in the Assembly resource from the same species, more than 15 Mbp, or is otherwise suspiciously long.
> * genome length too small - total non-gapped sequence length of the assembly is less than half that of the average for the genomes in the Assembly resource from the same species, less than 300 Kbp, or is otherwise suspiciously short.
> * hybrid - sequences from a hybrid between different species, strains or isolates.
> * low quality sequence - long stretches of the sequence have a high proportion of ambiguous bases, are low complexity, or have some other indication that the sequence quality is low.
> * mixed culture - sequences come from two or more organisms that were not cultured separately.
> * misassembled - alignment to related genome assemblies or other evidence indicates the assembly is likely to have errors.
> * partial - the assembly has only partial genome representation.
> * sequence duplications - assembly has one or more large duplications.
> * unverified source organism - the origin of the assembly is misidentified.

---
Tags: genome, data-download, ncbi
---


title: Parallelize or qsub a bash script

text:
# Question
Title: Parallelize or qsub a bash script
I have a bash script that I would like to parellelize to run on multiple nodes. My goal is to run my python sample\_script.py script on pairwise comparisons of samples to see if their variants are a match or not (output is written by echo for the function so we can compare results). I unfortunately can’t post the contents of the Python script, but it is basically comparing germline variants between samples to determine if a sample is a match or not. I’m trying to figure out how to make the script run faster or run multiple pair wise comparisons at once then concatenate results. The samples I'm comparing are listed in the array. Running the script on multiple samples takes too long, so I'd like to qsub the process or run the script on multiple nodes, then concatenate the results. Can someone advise on the best way to tackle this? I don't have much experience with jobs using qsub. Should I consider an SGE array job and if so how would I go about this?

Example code below using bash within jupyter notebook:

```
%%bash
array=( A037375001 A037646101 A037710601 )
for i in "${array[@]}"
      do
        for j in "${array[@]}"
        do
        python3 sample_script.py --run-sample-id-A "$i" --run-sample-id-B "$j"
        echo run-sample-id-A "$i" run-sample-id-B "$j"
      done
    done > "/outputdir/out.csv"

```

Thank you in advance!

# Answer
> 5 votes
What's best is quite subjective. My preference would be use Nextflow to abstract away the underlying job submission system and then use the sge executor to submit the jobs to the cluster. The only issue I see, is that it looks like you're feeding in a list of sample IDs, which your script presumably uses to find some VCF files containing the variants. If the script could be refactored to accept file inputs, we can have Nextflow localize them into its working directory. This would ensure your workflow (however trivial) is portable and could be run on AWS Batch, for example. The following example uses the combine operator to create the cartesian product and the collectFile operator to gather the results. With the following in a file called `main.nf`:

```
params.vcf_files = './path/to/*.vcf'

process compare_variants {

    input:
    tuple(
        path(sample_A_vcf, stageAs: 'A/*'),
        path(sample_B_vcf, stageAs: 'B/*'),
    )

    output:
    stdout

    """
    your_script.py \\
        --sample-vcf-A "${sample_A_vcf}" \\
        --sample-vcf-B "${sample_B_vcf}"
    """
}

workflow {

     samples = Channel.fromPath( params.vcf_files )
     product = samples.combine(samples)

     compare_variants( product )
         .collectFile( name: 'results.txt' )
         .view { "Saved results to: $it" }
}

```

And the following in a file called `nextflow.config`:

```
process {

  executor = 'sge'
}

```

And your Python script in a folder called 'bin' in the root directory of your project repository. Nextflow automatically adds this folder to the $PATH in the execution environment. So make it executable, for example using `chmod +x your_script.py`. Then, with Nextflow installed give it a test:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.7
Launching `main.nf` [chaotic_watson] DSL2 - revision: ac224d863b
executor >  local (484)
[bc/a39b27] process > compare_variants (482) [100%] 484 of 484 ✔
Saved results to: /path/to/work/tmp/3e/95330f32edec3f109b197bc836b2ca/results.txt

```

# Answer
> 2 votes
Consider using a workflow language such as Snakemake, Nextflow or Cromwell to handle distributing the jobs for you. It will take a little more time investment upfront to set them up, but it will save you many hours and headaches later on.

# Answer
> 2 votes
You may use asub (array job submitter) for this. Modify your script to `myscript.sh`:

```
array=( A037375001 A037646101 A037710601 A037710701 A037753001 A037844701 A038153001 A038172201 A038184401 A038272502 A038298602 A038533901 A038797201 A038944801 A038973401 A039034701 A039048101 A039089601 A039164101 A039303301 A039323801 A039538301 )
for i in "${array[@]}"; do
  for j in "${array[@]}"; do
    echo "python3 sample_script.py --run-sample-id-A $i --run-sample-id-B $j > out-$i-$j.csv"
  done
done

```

Note that the main difference is this script prints command lines, not executing them. You then run asub with

```
bash myscript.sh | asub -q myqueue -j myjobname

```

asub will submit a job array. If you don't like asub, you may submit with

```
bash myscript.sh | xargs -i echo "echo '{}' | qsub -q myqueue" | sh

```

This will submit separate jobs, not a job array. It will be slower and less convenient but you don't need to install anything or learn anything.

With either approach, you need to run the following to collate all the results:

```
cat out-*-*.csv > our.csv

```

---

On a side note, this approach can be adapted for GNU parallel to spawn multiple processes on the same machine:

```
bash myscript.sh | parallel -j 16

```

# Answer
> 0 votes
**Background** Just for everyone else `SGE` = Sun Grid Engine and is the `qsub` system mentioned in the title of the question.

### `qsub`:

If by parallelising you mean multithreading a single job across loads of cores then the first thing is how many cores per node are there?

Usually it is around 10 (but could be much more). That exact number is super important to know because if you request more cores than exist on a single node the job will never run. At best qsub will refuse the job, at worst it will forever be stuck in the queue. For multi-processing on `qsub` its ...

```
qsub -pe omp 8  myscript.sh

```

Here its requesting 8 cores and myscript.sh will contain the shell command to run the python script. If there's a node with 10-cores 2 are in use, it will then load the job to give its max. capacity.

Thus, `qsub` does not do **course grain-parallelisation**, that is `MPI`. Thus `qsub` will only parallelise across the cores in a single `node`. This is without question a limitation for cluster computing.

`qsub` is more than just submitting to a queue, you need to monitor the queuing system pre- and post- submission. Firstly, this is for the number of existing jobs across that system. Secondly, this is to see whats happening to your stuff.

`qstat` is the way to understand whats there:

```
qstat -f # whats what on the cluster
qstat -q

```

This lists all the queues.

```
qstat -u username # looks at what a given user is doing
qstat long-queue

```

Targets availability on a specific queue.

A complete list of `qstat` is here. Also `qusage -l` is useful.

---

### "Parallelising" loops

If you are submitting loads of jobs which are all working independently in parallel it better to use a different strategy on `qsub`. It is better to submit each job separately to the queue via a `qsub` loop. Thus there would be two scripts the first is the script in the question, the second is a submission loop: this would comprise a `qsub` argument within a loop that would submit each job sequentially. Thus the array might be better in the `qsub` submission loop.

Monitoring the queue takes place as per normal via `qstat`.

#### Rationale

The reason for this is the way the queue works `qsub` prioritises single jobs over multi-threaded parallelisation, thus you get your results much quicker.

---
Tags: python, bash, parallel
---


# Question
Title: Can you infer orientation from vcf file data?
I have a dante labs dataset that I would like to cross reference with SNPedia but I am unclear about orientation. I've read the info on SNPedia and other places about it but I just don't get it yet. I realise this is my own stupidity, the explanations feel sound but they seem to be missing some key info that I can't find anywhere.

Anyway without being too tl;dr about what confuses me, I wonder if my dante vcf file has the orientation problem or not because it shows what the reference genome base is next to the sample base. So for example:

```
1 36571920 rs96067 G A

```

Chr:1, Pos: 36571920, rsid:rs96067, reference base:G, sample base: A, (another field later in the line defines that both bases are the alt base so this is an (A;A) where the reference had (G;G)).

On SNPedia the genotypes mentioned for rs96067 are C;C, C;T and T;T. Since the reference base in my dante file is a G, which does not show up as a possible reference base on the SNPedia page, does this mean that the dante file is flipped with respect to this gene compared to SNPedia and so I should flip one or other if I want to compare them (flip the base in my file or the data from SNPedia that is)? So I guess flip G-\>C and A-\>T meaning my AA becomes TT?

BTW what confuses me about the SNPedia explanation about orientation is this:

> Unfortunately if this was instead a SNP where the two alleles are A or T the same flipping logic falls down. We don't (yet) have a way to know for sure if you should be flipped or not, since both forms of your flip rs1234(A;A) and rs1234(T;T) are possible.

Why is AA not the same as TT in this case? Surely this just reflects the direction that the bases were read from, so detection of an AA implicitly means that if we had read from the other direction we would have read TT? I appreciate this is wrong but I don't understand why.

@terdon: Re where I was quoting from, https://www.snpedia.com/index.php/Ambiguous\_flip (https://www.snpedia.com/index.php/Orientation also has some info). I've been taking 'flip' as their word for inversion. Re the 'orientation problem' I meant the ambiguity re orientation mentioned in those pages. However you mentioned in your comment that the convention is to always use the + strand. I read up about that and understand now that VCFs use that convention, which answers my question (thanks!). Dbsnp states that it always uses +ve. So like can be compared with like, as long as you use SNP definitions for the correct reference build that your individual SNP VCF has been generated from, and you are using dbsnp as your database. The so called 'flipping' can happen between reference versions/builds and could trip you if you don't use the correct one. Thats what I have gathered to this point, do let me know if its not right.

On SNPedia however many rsid's seem to be defined as minus orientation. They have 2 fields, orientation and stabilised orientation. To quote from the above linked pages:

> Orientation indicates the orientation reported in the most current build (currently GRCh38), which is reported below the genotypes in the Reference field. The chromosome number and nucleotide position in that reference build is shown next. StabilizedOrientation is the orientation that is relevant to the genotypes that have been defined in SNPedia for each SNP...

So they are saying that GRCh38 has both plus and minus orientation within and each genotype includes an orientation field? How does this fit with the standard '+ve' convention you mentioned? Maybe these are 2 different things and I just misunderstand. But it would seem far more sensible to simply adjust GRCh38 to all +ve. Is this 'orientation problem' just SNPedias problem? What exactly do they mean on the above linked pages on orientation and ambigous flips?

SNPedia seems pretty useful if it can be made to work - there is lots of associated data on each genotype page, the API is easy to use including for agregated queries etc, but it seems way less reliable than dbsnp despite apparently simply taking their data from dbsnp. Taking the above example I mentioned, rs96067 G\>A, dbsnp defines this well, for GRCh37 it shows NC\_000001.10:g.36571920G\>A. There is also a 'COL8A2 RefSeqGene NG\_016245.2:g.23766C\>T' entry which I don't fully understand. However on SNPedia as I said all that shows up are C;C, C;T and T;T genotypes. Maybe they are referring to the COL8A2 C\>A entry on dbsnp - or did they scrape this SNP def badly? Or are they showing transcript change rather than genomic as Ram RS suggested?

# Answer
SNPedia seems to show transcript base change and not genomic base change, so their C\>T is the same as the G\>A you note.

I don't understand SNPedia's allele notations completely, other sources are more straightforward. Does their `A;A` mean `A` on both strands? Then why do they talk about dbSNP showing `C;T` when dbSNP would in fact show `C>T` \- they seem to use `;` to mean alleles and base change interchangeably, and that is confusing. I recommend sticking to dbSNP for all information not found in SNPedia - something tells me SNPedia is not as trustworthy.
> 1 votes

# Answer
To answer (a small) part of my own question 'Why is AA not the same as TT in this case?' - I believe the answer to this is that the flanking sequence is important here. We can't just flip an AA to a TT in the middle of a strand without having any effect unless we also flip the flanking sequence also. What matters is the relation between the particular base in question and its flanking sequences. I didn't appreciate that before for some reason but its pretty obvious.
> 1 votes

# Answer
SNPedia introduces its own orientation to cut the tie to a reference genome, such that the bulk data is independent of the reference genome in use. The intention is good. The implementation is ok. It is just a bit confusing.

First of all, orientation is always relative; there is no absolute orientation. When you talk about plus or minus, you have to set a "reference" which could be a reference genome, transcription orientation or even something arbitrary.

Then VCF. Orientation in VCF is always relative to the reference that is used to generate the VCF.

In case of SNPedia, "Orientation" is secondary information that is relative to GRCh38, the current human reference genome (PS: this field changes with the reference genome but the internal primary strand system doesn't). rs96067 has two alleles C and T. Because its "Orientation" is "minus", the alleles on the GRCh38 strand is G and A, consistent with the dbSNP information. I checked a few other SNPs. It seems that SNPedia doesn't tell you which is the allele in the reference genome. This is actually fine because what matters is the genotype.

"Stabilized" orientation is relative to an arbitrary orientation stored in the SNPedia database. It is not directly relevant to your question when both you and SNPedia use the same reference genome.

How to use SNPedia? When you have a VCF in the GRCh38 coordinate, traverse rsIDs in the VCF. For each rsID, extract the genotype in the GRCh38 coordinate. Then get the SNPedia "Orientation" and "Geno" based on the rsID. If the orientation is "minus", complement "Geno" (i.e. "A;C" to "T;G", etc). Compare the genoptype from the VCF to the genotype association from SNPedia to predict the effect of genotype.

---

PS (clarification in response to new comments): if you design a database that ties to GRCh38 from ground up, you may need to flip many genotypes when the reference genome is changed. This process is complex and error prone. With their current design, if the reference genome is changed, they only need to update "Orientation" without changing the core data. This is a thoughtful design, though the UI could be made clearer.
> 1 votes

# Answer
Strand flips are common when using different VCF generation software. You've got one of the easier-to-determine variants there, because it's pretty clear that the orientation is incorrect between SNPedia and your VCF file. Illumina is particularly frustrating in this regard because they have a "top"/"bottom" definition, which can sometimes be inconsistent with the "reference positive" convention.

The reason strand for complementary variants can't be \[easily\] determined is because the usual way to fix strand in a bulk fashion is to check for which of the variant or it's complement are present in the reference database. If the variant is a complementary variant, then *both* will be found (one as the reference allele, and the other as the alternate variant allele). There are some statistical ways to work round this (e.g. looking at the most common variant in your dataset, and comparing it to the highest allele frequency in the reference database), but those methods will fail if the ancestral background of your study population is substantially different from the reference population (which was the case for the population genetics work I was doing).

The ultimate reference for RS numbers is dbSNP, because that's where the RS numbers are dished out for new SNPs.

Illumina claims that their top/bottom definition allows them to precisely determine complementary strand orientation, but it's useless unless the full sequence context is provided (i.e. the surrounding sequence up until strand can be determined by reverse complementation). Here's the explanation from Illumina:

> For SNPs that are not \[A/T\] or \[G/C\], A is always on the top strand and T is always the bottom strand. A and T nucleotides are the “A alleles”; G and C nucleotides are the “B alleles”. If the SNP is \[A/T\] or \[G/C\]: Use sequence walking to determine TOP/BOT strands, then assign A/B alleles.

```
rs903997 ...ATGAGA*A*AGT[C/G]TGA*G*AGTGCA...

```

By Illumina's definition, the above sequence is designated "top", because if you walk outwards from the variant location, you eventually get to a situation where A is present on one side, and C or G is present on the other side. The reverse complement of this sequence would be designated "bottom":

```
rs903997 ...TGCACT*C*TCA[C/G]ACT*T*TCTCAT...

```

This is explained more exhaustively in Illumina's technical note:

https://www.illumina.com/documents/products/technotes/technote\_topbot.pdf

When I was puzzling over this, I gave up trying to work this out, and filtered out all complementary variants in my dataset.
> 0 votes

---
Tags: snp
---


# Question
Title: How to extract sequences of a specific gene from genomes, with genomes' .fna, .gff and results from KofamScan?
I am trying to do alignment among orthologs and paralogs of a gene but I encountered some problems. I have came up with a workflow, but the sequences I extracted are not correct. They are of hugely varying length, and the alignment is very poor.

I have attached the workflow below, could anyone suggest if I am missing out on anything or if I am interpreting the files wrongly? Really appreciate your help.

My query is: with a list of target genes and bacteria genomes (.fna), for each gene, extract all sequences of it from all genomes and put together into one file for alignment

# Answer
> 1 votes
The answer is very simple, just align from amino acids. Amino acids are not subject to saturation - which is what you are encountering and is smashing up the nucleotide alignment. Amino acid will result in a high quality alignment. Thus you think the sequences are wrong, I think the sequences are simply divergent.

The other issues of variation in length - that depends on the how you've assembled the genes and exactly what project this is. *Sinorhizobium* is a pretty specific bacteria, so it's pretty clear what the biology is - its quite possibly legume agriculture. However, the precise molecular study is not so clear. Basically, what I'm saying is 'variations' in length are not necessarily 'errors' if this is metagenomics.

You could extract the amino acids from your `fna` file. Or ... use gff format to get the amino acids via `gffread -y protein.fa` that sort of thing.

The following might not make much sense ... but hope it helps ... to obtain the DNA alignment ... you convert the amino acid alignment back to a nucleotide alignment by transferring the alignment from amino acids on to the nucleotide data using Emboss' `tranalign`, amongst many other programs.

The tranalign command is:

```
tranalign -asequence [nucleotide filename (fasta)] -bsequence [aligned protein filename (fasta)] -outseq

```

tranalign is also available online, albeit the program is nothing special and its just as easy to write the code.

You need a bit of code to translate the nucleotides to amino acids \[I use Biopython `translate()`\] or just use `seqotron` of `MEGAX` translate function (if you've a Mac) for a GUI.

**Warning** If you use this approach the taxa names (strain names) KOs etc ... need to be identical between amino acid and nucleotide files AND the order needs to be identical too. This can be a headache - but if you encounter errors here simply post back as a separate question and I'll forward the code to get around this 'glitch'.

I hope this all sort of makes sense.

---

HOWEVER, if you want to assess whether KOfamscan is giving the correct results (which is what you are questioning) you could use

KofamKOALA: KEGG Ortholog assignment based on profile HMM and adaptive score threshold available

This is a webserver available here that will do the same thing and should give the same results. The authors consider it an upgrade, but I've not used it.

---

There are ALOT of other approaches to this problem BTW and all have over-reliance on KEGG, whilst it is the most comprehensive database available, is not perfect. KEGG *per se* omits upwards of half of all genes. Nobody else needs to know that but as long as you are aware of this thats fine.

---
Tags: gff3, kegg
---


# Question
Title: Where can I find examples of published bioinformatics pipeline packages?
I am a software engineer and I am preparing a presentation to aspiring bioinformatics PhDs on how to use best-practice software engineering when publishing code (such as include documentation, modular design, include tests, ...). In particular my presentation will be focused on publishing "pipeline" code, that is, code that is mainly focused on transforming data to a suitable shape for analysis.

I am trying to find good example of published bioinformatics pipelines that I can point students to, but as I am not a bioinformatician I am struggling to do so, so would like your help. It doesn't matter if the published pipeline is super-niche or not very popular so long as you think it is engineered well so that I can look it up.

Specifically the published code should have: adequate documentation, testing methodology, modular design and be easy to install.

# Answer
> 4 votes
I think the nextflow core (https://nf-co.re) pipelines should satisfy most of your requirements. These pipelines are written in nextflow (https://nextflow.io) which allows scaling from personal to hpc infrastructure rather efficiently. The core pipelines are usually quite well documented, and have a module design. I have not really looked into the testing methodology, but I think the modules themselves are tested. Running them is easy once nextflow is installed on your system.

# Answer
> 4 votes
The GATK team maintains a GitHub repository with their pipelines: https://github.com/gatk-workflows. They use the Cromwell Workflow Management System, an open source tool created by the Broad Institute.

You are probably most interested in https://github.com/gatk-workflows/gatk4-data-processing, their workflow for aliging fastq files against a reference genome to produce a bam file which can be used for downstream analyses such as variant calling:

> # gatk4-data-processing
>
> ## Purpose :
>
> Workflows for processing high-throughput sequencing data for variant discovery with GATK4 and related tools.
>
> ## processing-for-variant-discovery-gatk4 :
>
> The processing-for-variant-discovery-gatk4 WDL pipeline implements data pre-processing according to the GATK Best Practices. The workflow takes as input an unmapped BAM list file (text file containing paths to unmapped bam files) to perform preprocessing tasks such as mapping, marking duplicates, and base recalibration. It produces a single BAM file and its index suitable for variant discovery analysis using tools such as Haplotypecaller.
>
> * If you are starting with FASTQ files visit the seq-format-conversion repository for workflows to convert FASTQs to unmapped BAMS.
> * The processing-for-variant-discovery-gatk4 provides quick and general processing for sequence data using the latest releases of GATK. If users are interested in a more elaborate version of this workflow with quality control tasks and routinely tested for validity (useful in production environments) then visit the gatk4-genome-processing-pipeline repository.
> * The BAM output from processing-for-variant-discovery-gatk4 can be used to perform a variety of other analysis like somatic short variant discovery, germline short variant discovery, or germline copy number variant discovery. Visit the GATK Best Practices documentation to determine what to do next with the BAM files.

---
Tags: data-preprocessing, data-retrieval, publishing
---


# Question
Title: How to differentiate DNA fastq and RNA fastq files?
I have 2 sets of fastq files from another collaborator. The first is contains exome data and the second contains RNAseq data. But both have the RNA in the name but a different ID. How do I differentiate between a fastq file between exome seq and RNAseq?

# Answer
welcome to the Bioinformatics StackExchange! This is one of those cases where the best solution is to just ask your collaborators. Don't bother with anything technical from the data itself because it's only going to waste your time.
> 5 votes

# Answer
Absolutely agree with James that you should ask your collaborators, but purely as an academic exercise... assuming you have a good reference genome+annotation (e.g. human/mouse), whole exome seq should align to mostly annotated exonic, RNA-seq to exonic+intronic, and whole genome seq to exonic+intronic+intergenic sequences. Sequence duplication levels may also be higher for RNA-Seq with highly expressed genes often 'deeply' sequenced
> 3 votes

---
Tags: rna-seq, fastq, exome
---


# Question
Title: GATK g CNV code Aborts Before Running Main Code
I have been struggling to get the code I wrote for running GATK g CNV on my Linux machine (Steps Here: https://gatk.broadinstitute.org/hc/en-us/articles/360035531152 ) to work. Specifically it seems to abort a run before running the main code (I even removed my DetermineGermlineContigPloidy step to no avail). Below is my code and then the error:

```

## Paths To Be Set

echo "My paths have been set"

### Set Bam Directory:
#### Make sure to set a seperate directory for the batch you are doing

BAM_DIR=/pathToThe/Batch_I_Am_Running

### Set the path to the reference genome:
reference=/pathToMy/primary_assembly.fa

#model_dir =/path/to/my/GATK-gCNV

### Set the path to the singularity image:
SINGULARITY_CACHEDIR=/path/to/my/GATK
singularity_image=$SINGULARITY_CACHEDIR/gatk_latest.sif

### Set Number Of Bins
#### Normally 100 to 200*

NUM_BINS=200

### Setting the number of threads ?? Like Cores ??

NUM_THREADS=20

#______________________________________________________________________________________________________________________________________________________________

## Main code

echo "I am about to start the main code"

for c in "${CHROMOSOME_NAMES[@]}"; do

  # loop over the directories containing the BAM files
  for dir in "${BAM_DIR}"/*/; do

    # get the name of the directory without the path
    dir_name="$(basename "${dir}")"

    # get the name of the BAM file without the extension
    bam_name="$(basename "${dir}"/*.bam .bam)"

    echo "Processing file: $bam_name"

    echo "I am starting Step 1: Preprocessing"

    # Step 1: Preprocessing
    singularity exec $singularity_image PreprocessIntervals \
  -R ${reference} \
      -I ${dir}/${bam_name}.bam \
      -L targets.interval_list \
      --bin-length 0 \
      -imr OVERLAPPING_ONLY \
      -O targets.preprocessed.interval_list

    echo "I am starting Step 2: AnnotateIntervals with GC content"

    # Step 2: AnnotateIntervals with GC content
    singularity exec $singularity_image AnnotateIntervals \
  -L chr${c}.interval_list \
      -I ${dir}/${bam_name}.bam \
      -R ${reference} \
  -imr OVERLAPPING_ONLY \
  -O chr${c}.annotated.tsv

    echo "I am starting Step 3: CollectReadCounts"

    # Step 3: Collect read counts for the samples
    singularity exec $singularity_image CollectReadCounts \
  -I ${dir}/${bam_name}.bam \
  -L targets.preprocessed.interval_list \
  --interval-merging-rule OVERLAPPING_ONLY \
  --format TSV \
  -O ${dir_name}.${bam_name}.counts.tsv

    echo "I am starting Step 4: FilteringIntervals based on GC-content & read depth"

    # Step 4: FilterIntervals based on GC-content and read depth
    singularity exec $singularity_image FilterIntervals \
  -L chr${c}.interval_list \
      --annotated-intervals chr${c}.annotated.tsv \
  -imr OVERLAPPING_ONLY \
  "${dir}"/*.counts.tsv \
      --max-panel-norm 1.5 \
      --min-panel-norm 0.3 \
      -O chr${c}.filtered.interval_list

    echo "Not doing Step 5 DetermineGermlineContigPloidy for now"

    # Step 5: DetermineGermlineContigPloidy (Mandatory Step)
   # singularity exec $singularity_image DetermineGermlineContigPloidy \
 # --model-parameters ${model_dir}/model_parameters.txt \
     # --contig-ploidy-calls ${model_dir}/contig_ploidy_calls.tsv \
 # --contig-ploidy-priors chr20XY_contig_ploidy_priors.tsv \
 # --output-prefix ${model_dir}/ploidy \
     # --standardized-output \
     # $(find "${output_dir}" -name "*.calls.tsv" -printf '%p ') \

#    for file in ${input_dir}/*.tsv
#    do

    # Get sample ID from filename
#      sample=$(basename ${file} .tsv)
#      # Run DetermineGermlineContigPloidy command
#      singularity exec ${singularity_image} DetermineGermlineContigPloidy \
#        --input ${file} \
#        --contig-ploidy-priors ${contig_ploidy_priors} \
#        --output-prefix ${output_dir}/${sample}_ploidy \
#        --standardized-output
#    done

    echo "Creating shards via bins"

    # Creating shards via bins
    singularity exec $singularity_image IntervalListTools \
  SplitIntervals \
  -R ${reference} \
      -L chr${c}.interval_list \
  -scatter ${NUM_BINS} \
      -o ${dir}/${c}_shard.interval_list

    echo "Starting the Shard Loop that will run GermlineCNVCaller and PostprocessGermlineCNVCalls"

    # Loop over the shards and run GermlineCNVCaller and PostprocessGermlineCNVCalls
    for shard_file in "${dir}/${c}"*_shard.interval_list; do

      echo "Starting Step 6: Run GermlineCNVCaller (Mandatory Step)"

      # Step 6: Run GermlineCNVCaller (Mandatory Step)
      singularity exec $singularity_image GermlineCNVCaller \
    --run-mode COHORT \
    -L "${shard_file}" \
        -I $(find "${dir}" -name "*.tsv" -printf '%p ') \
        --minimum-interval-median-percentile 5.0 \
        --maximum-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --output-prefix $(basename "${shard_file}" .interval_list)

      echo "Starting Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)"

      # Step 7: Run PostprocessGermlineCNVCalls (Mandatory Step)
      singularity exec $singularity_image PostprocessGermlineCNVCalls \
    --model-shard-path "${shard_file}" \
        --cnv-call-path $(basename "${shard_file}" .interval_list).calls.tsv \
        --output-dir postprocess-sm \
        --scatter-count ${NUM_THREADS} \
        --num-scatter-plot-columns 8 \
        --max-zeros-in-sample-percentage 40.0 \
        --annotation-cutoffs "QD=2.0" "SOR=2.0" "FS=10.0" "MQ=50.0" \
        --interval-merging-rule OVERLAPPING_ONLY \
        --minimum-contig-length 46709983 \
        --minimum-interval-median-percentile 5.0

    done

  done

done

    #______________________________________________________________________________________________________________________________________________________________

```

The ".out" file I received after the run aborted is:

> My paths have been set
>
> I am about to start the main code

Why would it abort so soon (I changed the number of Gs used and that did not have any effect)?

# Answer
You are using `${CHROMOSOME_NAMES[@]}` but you never define that, so when your code gets here:

```
for c in "${CHROMOSOME_NAMES[@]}"; do
...
done

```

It doesn't enter the `for` loop because `"${CHROMOSOME_NAMES[@]}"` is an empty array. It's the same thing as doing this in a shell session:

```
$ for c in "${CHROMOSOME_NAMES[@]}"; do echo "This is chromosome $c"; done
$

```

Just like your script, that will produce no output since there is nothing for it to loop over.
> 1 votes

---
Tags: gatk, cnv
---


# Question
Title: How to calculate the Maximum Ladder Distance (MLD) of ssRNA from Dot-Bracket notation?
The ladder distance (LD) between residue i and residue j is the minimum number of base pairs (created due to secondary structure of the ssRNA) one has to cross when going from i to j.

The Maximum Ladder Distance of a ssRNA molecule is the maximum ladder distance between any two of the residues of this molecule\[\*\].

Given the secondary structure of a ssRNA molecule in dot-bracket notation (e.g. by RNAfold), **how can I calculate the Maximum Ladder Distance of this molecule**?

Below is an illustration from \[\*\]. The MLD is 11, since this is the number of basepairs on the yellow line.

Edit: The dot-bracket notation for this molecule would be

```
..(((((.......))...((.((((((.......))))))..)).)).

```

\*: Yoffe AM, Prinsen P, Gopal A, Knobler CM, Gelbart WM, et al. (2008) Predicting the sizes of large RNA molecules. Proc Natl Acad Sci U S A 105: 16153–16158

# Answer
There are different implementations for calcuating the $\rm MLD$. One possible approach is to consider $\rm MLD$ as the *heighest value of all possible mountain plots*, generated by shifting the beginning of your sequence. The highest value of a mountain plot tells you the maximal ladder distance that can be reached from *the first nucleotide*. But the longest path might not begin at the first nucleotide, and that is why you have to go through all possible shifts of the sequence.

Here is an example of this implementation in Python. The algorithm is as following:

1. Convert dot-bracket notation to a list of base pairs (`dot2pairs`).
2. Calculate the first mountain plot.
3. Calculate required `shifts`. *This is not strictly neccessary because you can swipe through all possible `shifts = list(range(0, length))`, but it gives better performance.*
4. Iterate through all `shifts`, calculate the `MLD` at each step, and return the maximal one.

```
import numpy as np
from collections import deque

def dot2pairs(dot):
    stack, pairs = deque(), []

    for i, n in enumerate(dot):
        if n == '(':
            stack.append(i)
        elif n == ')':
            pairs.append((stack.pop(), i))

    return np.array(pairs)

def dot2MLD(dot):
    # Convert to a list of base pairs
    pairs0 = dot2pairs(dot)

    if len(pairs0) == 0:
        return 0

    MLD = 0
    shifts = [0]
    length = len(dot)
    seq = np.zeros(length, dtype=int)
    curr_mount = None

    while shifts:
        # Shift the sequence
        pairs = (pairs0 - shifts.pop()) % length

        # Find positions of opening and closing nucleotides
        opens, closes = np.min(pairs, axis=1), np.max(pairs, axis=1)

        seq.fill(0)
        seq[opens] = 1
        seq[closes] = -1

        mountain = np.cumsum(seq)

        MLD = max(MLD, int(np.max(mountain)))

        # In first iteration, calculate all required shifts
        if curr_mount is None:
            curr_mount = mountain[0]
            for i, m in enumerate(mountain):
                if m == curr_mount:
                    shifts.append(i)
                else:
                    curr_mount = m

    return MLD

```

Let's test it on your structure:

```
dot2MLD("..(((((.......)))...((.((((((.......))))))..)).)).")

```

> 11

And on a larger structure:

```
dot2MLD(s)

```

> 1141

It took ~ 5 seconds for this 17481 nt long structure (and it takes ~ 20 seconds if you don't precalculate requires `shifts`).
> 5 votes

---
Tags: rna
---


# Question
Title: Test for differences between groups of samples
Sorry if the answer to this should be obvious.

I have RNA-expression results from 24 samples which can be divided into 6 groups, (wildtype and two different mutants at two different ages) with a total of 24 samples. I have made a PCA plot and when describing it in my thesis said separation between samples is minimal. One of my examiners has said

> You state that the separation between samples in the principal component analysis is minimal, this needs to be quantified to differentiate between minimal but statistically significant, and not significantly different.

I have no idea how he expects me to test for significance. I know how to test for significant differences between individual samples, and for differences per gene/principle component. I don't know how to group samples together while testing for significance against multiple things or if I should use genes or principle components (probably the latter).

I'm writing it in R so if anyone knows which R function can do the kinds of statistical tests I need that would be a huge help. Alternatively is this something which I shouldn't really do and I should push back on?

I'm already performing significance testing using DESeq2. To the best of my knowledge that only tests for whether a single gene is significant. It doesn't test for whether the samples are significantly different as a whole which is what I'm after.

# Answer
PCA does not do any testing and has no p-values. You may state that a separation of less than `m` units in PC \> `i` is not significant while a separation of greater than `n` units in PC \< `i` is significant as PCs \< `i` explain X% of the variation (`i` being 2 or 3). That's how you can back your statement on the separation being minimal.

That is, in this sample plot:

PC2 distance of \< 0.05 might be called trivial compared to PC1 distance \> 0.5 if PC1 explains 80% of the variance. This is just 2 PCs though, so your actual case will need to be more robust.
> 2 votes

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


# Question
Title: How to get enriched pathways in the data using continous statistic measure?
I was doing `pathway enrichment analysis` using the below code

```
p_value = scipy.stats.hypergeom.sf(len(intersection) - 1, geneset_size, len(pathway), len(hits))
        if p_value < 0.05:
            for gene in list(pathway):
                results = results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

```

where,

```
hits = df["Gene"][np.logical_and(df["Model"] == model, df["p_value_anova"] < 0.05)]
intersection = pathway.intersection(hits)
Model = model name (like A, B, C) in input file
Gene = Gene from the input gene
Enriched pathway=Pathways that are enriched from the Pathway input file provided

```

Now the problem is that earlier this checking pathway enrichment on the basis of p\_value in `hits` was working but now using cutoff (0.05) is resulting in loss of all pathways. I am looking to test the enrichment using the statistical test that will work on continous measure than binarizing the P-value.

# Answer
> 1 votes
Just to mention I personally think there is a higher loop above the loop stated above.

I've not looked at each bit of code in detail but there's a bug is on line 4:

```
results = results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

```

Firstly, this is fashionable coding because you are dumping straight into a `pandas` dataframe (described below). The correction is simply ... declare `pandas` in advance

```
import pandas as pd
results = pd.DataFrame()
p_value = # code missing
if p_value < 0.05:
    for gene in list(pathway):
        results.append({"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}, ignore_index=True)

```

The thing is if there is a higher loop above `list(pathway)` then `results = results.append(....` will over-ride the last set of results. This would explain why results 'disappear' because if the last data in the loop is a blank - or very little output, all previous stuff gets deleted.

I personally would use a dictionary here. Thus

```
# declare the dictionary before all loops
results = {"Gene": [], "Model": [], "Enriched Pathway": []}
p_value = # code missing
if p_value < 0.05:
    for gene in list(pathway):
       genedict = {"Gene": gene, "Model": model, "Enriched Pathway": Pathway_name}
       for k,v in mydicti.values():
           results[v].append(v)

```

... you could then transfer this to `pandas` once the loops are complete via,

```
df = pd.DataFrame(genedict)

```

I do suspect there's other issues, but as a general rule `results = results.append(...` is better declaring `results` first and then appending to `results`, unless you've a very specific reason (e.g. writing each loop into a separate file). The `pandas` thing is kinda clever for big data sets, because `pandas` is fast, whilst dictionaries are slow - for big data (like \>\>500 000 big). The interchangeability between `dictionaries` and `pandas` means it's kinda easier conceptually (and code-wise) to switch at the coders convenience in my personal opinion.

---
Tags: python, statistics, pathway, go-enrichment, pandas
---


# Question
Title: How can I export all and only the publications cited by a specific article?
I want to do a systematic literature review of genetic variant annotation. As part of this I am interested in all the databases available to annotate genetic variants.

I want to export all the references cited by the article with PubMed ID 35943817 (https://pubmed.ncbi.nlm.nih.gov/35943817/), because many of these references are relevant to my review. I can see that they number 93. I want them in a CSV table with their PubMed ID, Title, etc., similar to CSV tables obtained in a search after clicking "Save citations to file", Format "CSV".

# Answer
Something like this will work, it did just now for me. You might want to extract different info from the records though.
This exports to a table (bibliography.tsv) columns containing:

* pmid
* authors
* title
* journal
* year

```
max_attempts = 5
def fetch_pubmed_data(pmid):
    Entrez.email = "you@email.com"  # Replace with your email address
    handle = Entrez.efetch(db="pubmed", id=pmid, rettype="xml", retmode="text")
    record = Entrez.read(handle)
    time.sleep(1)
    return record

def extract_cited_articles(pmid):

    record = fetch_pubmed_data(pmid)
    cited_pmids = []

    record_article = record['PubmedArticle']

    if "PubmedData" in record_article[0] and "ReferenceList" in record_article[0]["PubmedData"]:
        for reference in record_article[0]["PubmedData"]["ReferenceList"][0]["Reference"]:
            if "ArticleIdList" in reference and len(reference["ArticleIdList"]) > 0:
                cited_pmids.append(reference["ArticleIdList"][0])

    return cited_pmids

def extract_bibliographic_info(pmid_list, max_attempts):
    bibliographic_info = {}
    bibliographic_info_ls = []

    for pmid in pmid_list:

        authors = title = journal = year = ""
        record = {}

        for attempt in range(max_attempts):
            try:
                record = fetch_pubmed_data(pmid)
                break

            except:
                time.sleep(5)

        if 'PubmedArticle' in record:
            if len(record['PubmedArticle']) >0:
                article = record['PubmedArticle'][0]["MedlineCitation"]["Article"]

                if "AuthorList" in article:
                    authors = ", ".join([f"{author['LastName']} {author['Initials']}" for author in article["AuthorList"] if 'LastName' in author and 'Initials' in author])
                    title = article["ArticleTitle"]
                    journal = article["Journal"]["Title"]
                    year = article["Journal"]["JournalIssue"]["PubDate"]["Year"]

        print (pmid, authors, title, journal, year)
        bibliographic_info_ls.append([pmid, authors, title, journal, year])

    return bibliographic_info_ls

def main():
    pmid = input("Enter the PubMed article ID: ")
    cited_pmids = extract_cited_articles(pmid)
    bibliographic_info_ls = extract_bibliographic_info(cited_pmids, max_attempts)

    df = pd.DataFrame(bibliographic_info_ls, columns=["pmid", "authors", "title", "journal", "year"])
    df.to_csv("bibliography.tsv", sep="\t", index=False)
    return df

if __name__ == "__main__":
    df = main()

```
> 2 votes

---
Tags: literature-search
---


# Question
Title: Combine separate fasta files from a list.txt
I have about 10 fasta files that I need to pull from a giant folder of other fasta files (has about 100 files). I have a list of the files I need in a .txt file separated by line and I want to pull out the 10 fasta files listed and combine them into one big file.

I tried this command at first: `xargs < names.txt cat > names.fasta`

But it gave me this error: `cat: 'read_1.fasta'$'\r': No such file or directory`

For all the files except the last read in the list which it managed to add to the `names.fasta` .

I'm not sure what the error message means given I can find the files in my folder so they should be there. I can do this manually but I know it'll take a while to do.

**EDIT**: I changed the list to be separated by space and it actually pulled out all the files I needed! However, is there a way to account for reading in the names by newline instead of space??

# Answer
> 2 votes
This should work

```
xargs < names.txt cat > names.fasta

```

Or else try,

```
cat names.txt | xargs cat > newdir/names.fasta

```

What I think is within `names.txt` there's a typo or a missing file and that's what it's complaining about.

Perhaps try using

```
ls > tmp
cat names.txt | sort > tmp2
diff tmp tmp2

```

If there are any differences that should pick it up.

# Answer
> 1 votes
Just to address an edge case, run `dos2unix names.txt` first and then a modified version of @M\__'s command

```
cat names.txt | xargs -I v_file cat v_file >> names.fasta

```

Based on your edit, it looks like your working platform has different line endings compared to the platform on which `names.txt` was last edited/saved.

---
Tags: fasta, command-line
---


# Question
Title: how to handle variable output from nextflow
I have this output which is going to do bwa next which is from join operator

```
[DF122_SO11378_A03, DF122_SO11378, DF122_SO11378_A03.norm.merged.fq, DF122_SO11378_A03.final.gfa, DF122_SO11378_A03.final.fasta]

```

in this case, `pair_id` is `DF122_SO11378_A03`, `plate_id` is `DF122_SO11378`, `fq` is `DF122_SO11378_A03.norm.merged.fq`, `DF122_SO11378_A03.final.gfa` is `gfa` `DF122_SO11378_A03.final.fasta` is `fa`

I have following code like this

```
process bwa {

     input:
     tuple val(pair_id), val(plate_id), path(fq), path(gfa), path(fa) from align_in

     output:
     path "*bam" into bam
     tuple val(pair_id), val(plate_id), path ("*.csv") into metrics

     """
     if [ -s $fa ]; then
     bwa index $fa

     bwa mem $fa $fq \
      | samtools view -bh -F2048 - \
      | samtools sort > ${pair_id}.bam

     samtools depth -a ${pair_id}.bam > ${pair_id}.depth.csv
     samtools view -c ${pair_id}.bam >${pair_id}.count.csv

     else
     touch ${pair_id}.depth.csv
     touch ${pair_id}.count.csv
     touch ${pair_id}.bam
     fi
     """
}

```

However, now my upstream work flow changed a little, for some sample they have more than one gfa and fasta output like this

```
[DF122_SO11378_A04, DF122_SO11378, DF122_SO11378_A04.norm.merged.fq, [DF122_SO11378_A04.2.final.gfa, DF122_SO11378_A04.final.gfa], [DF122_SO11378_A04.2.final.fasta, DF122_SO11378_A04.final.fasta]]

```

in this case, `pair_id` is `DF122_SO11378_A04`, `plate_id` is `DF122_SO11378`, `fq` is `DF122_SO11378_A04.norm.merged.fq`, `[DF122_SO11378_A04.2.final.gfa, DF122_SO11378_A04.final.gfa]` is `gfa`, `[DF122_SO11378_A04.2.final.fasta, DF122_SO11378_A04.final.fasta]` is `fa`

But some samples still have one gfa and fasta as shown earlier.

How can I adapt my bwa process to handle such situations?

I am using DSL1.

# Answer
If you're looking for the transposition, you could use the `transpose` operator to get what you want. But to handle the variable output, you would need to first check what types you have. For this, you could use something like the following for example:

```
input_ch
    .map { pair_id, plate_id, fastq, gfa, fasta ->

        def gfa_files = gfa instanceof Path ? [gfa] : gfa
        def fasta_files = fasta instanceof Path ? [fasta] : fasta

        tuple( pair_id, plate_id, fastq, gfa_files, fasta_files )
    }
    .transpose()
    .view()

```

Note that the old DSL1 is no longer supported.
> 1 votes

---
Tags: nextflow
---


# Question
Title: Singularity & GATK 4.4.0.0. CreateSequenceDictionary
I ran the following code in bash (I added print statements for ease):

```
# Module loading

echo "I am loading the specified modules"

module load Singularity/version_number

# Define the path to the singularity:

singularity_image="/path/to/folder/GATK/gatk_latest.sif"

echo "I am going to make the `.dict ` file"

singularity exec $singularity_image /gatk/gatk --java-options -Xms8000m \
  CreateSequenceDictionary \
  -R /path/to/my/primary_assembly.fa \
  -O /path/to/the/primary_assembly.dict

echo "I am all done"

```

However, I get the following `out` file:

```
I am loading the specified modules
job_scripts_number: line 23:
.dict: command not found I am going to make the  file
INFO  NativeLibraryLoader - Loading libgkl_compression.so from
jar:file:/gatk/gatk-package-4.4.0.0-local.jar!/path/to/libgkl_compression.so [Today's date] CreateSequenceDictionary --OUTPUT
path/to/my/primary_assembly.dict
--REFERENCE /path/to/my/primary_assembly.fa
--TRUNCATE_NAMES_AT_WHITESPACE true --NUM_SEQUENCES 2147483647 --VERBOSITY INFO --QUIET false --VALIDATION_STRINGENCY STRICT --COMPRESSION_LEVEL 2 --MAX_RECORDS_IN_RAM 999999 --CREATE_INDEX false --CREATE_MD5_FILE false --help false --version false --showHidden false --USE_JDK_DEFLATER false --USE_JDK_INFLATER false [Today's date] Executing as user@node_number on Linux version
; OpenJDK 64-Bit Server VM version;
Deflater: A; Inflater: A; Provider is available; Picard
version: Version:4.4.0.0 [Today's date]
picard.sam.CreateSequenceDictionary done. Elapsed time: 0.53 minutes.
Runtime.totalMemory()=8388608000 Tool returned: 0 Using GATK jar
/gatk/gatk-package-4.4.0.0-local.jar Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -Xms8000m -jar /gatk/gatk-package-4.4.0.0-local.jar CreateSequenceDictionary -R
/path/to/primary_assembly.fa
-O /path/to/primary_assembly.dict
I am all done

```

Why would I get the error ".dict: command not found I am going to make the file " especially when I look at the `head` of the .dict file and it seems to look the way it should?

# Answer
This has nothing to do with Singularity/Docker or GATK. You get the error `.dict: command not found` because you've used backticks to perform a command substitution in a prior `echo` command:

> `echo "I am going to make the `.dict ` file"`

Backticks are legacy syntax, you could have also written:

```
echo "I am going to make the $(.dict) file"

```

But neither of these makes any sense - there is no `.dict` command. Of course single quotes would prevent the command substitution:

```
echo 'I am going to make the `.dict` file'

```

But this might cause some confusion when the output is displayed. A better choice would be to use double quotes, but without the backticks. For example:

```
echo "I am going to make the '.dict' file"

```
> 5 votes

---
Tags: gatk, docker, dictionary, singularity
---


# Question
Title: In CellRanger output, what are the reads with an xf:i:17 tag?
When running CellRanger on 10x scRNA-Seq data, a BAM file is output with the aligned reads. These reads have an `xf` tag, described here:

> The bits of this tag are interpreted as follows:
>
> * 1 The read is confidently mapped to the transcriptome
> * 2 This read's barcode, UMI, and feature combination was discarded in favor of a different feature with higher read support
> * 4 This read pair maps to a discordant pair of genes, and is not treated as a UMI count
> * 8 This read is representative for a molecule and is treated as a UMI count
> * 16 This read maps to exactly one feature, and is identical to bit 1 for transcriptomic reads. Notably, this bit is set for a Feature Barcode read, while bit 1 is not
> * 32 This read was removed by targeted UMI filtering.

In a real dataset, I have this distribution of xf values:

```
$ samtools view possorted_genome_bam.bam |  egrep -o "xf:i:[0-9]+" | sort | uniq -c
66263599 xf:i:0
99892363 xf:i:17
  245179 xf:i:19
 8283967 xf:i:25

```

It is clear to me that `0` are unaligned, `19 = 16+2+1` are discarded, `25 = 16+8+1` are kept (and used in the feature\_barcode\_matrix); but I do not understand what happens to the `17 = 16 +1`. They seem to be good reads, they are aligned on genes, yet they are discarded.

One possible explanation from my reading of this paragraph is that the reads with `17` are PCR duplicates, thus ignored, whereas the reads with `19` are in a group matching different genes (several reads with the same BC/UMI but matching different genes). But this doesn't look right when exploring the BAM file manually.

# Answer
> 1 votes
From the documentation at https://kb.10xgenomics.com/hc/en-us/articles/115003646912-How-is-sequencing-saturation-calculated-

> All of these samflag 0x400 reads have an xf tag value of 17, which consist of bits 1 and 16. This also means these reads do not have the xf bit of 8, which mark representative reads from a group of duplicates. The converse isn't true though. The xf17 consist mostly of samflag 0x400 duplicate reads but also of samflag nonduplicate reads.

---
Tags: sam, 10x-genomics, cellranger
---


# Question
Title: What does it mean by a Local Coordinate in the case of protein chains?
> In Figure 3a, each residue has its own local coordinates, which are illustrated shown using different colors.

Here are two more papers:

---

By global coordinate, I understand the coordinate system in which we view the entire protein chain.

What does it **generally** (not only in this particular paper) mean by a Local Coordinate in the case of protein chains?

Why was it necessary to introduce this terminology when they could have just used the term "vector"?

Can you supply me with an example with an explanatory diagram?

# Answer
> 1 votes
For proteins and molecules in general, and within the meaning of this paper, "internal coordinates" are any representation that are not dependent on the molecule rotation or translation in a Cartesian space.

Calling them "internal coordinates" may be slightly confusing, and a more correct term could be "internal representation"

That means, I guess, that internal coordinates **could** be, as you mention, a list of vector in a Cartesian space, where the Cartesian space "orientation" would be decided by, say, principle component analysis (or any other arbitrary rule depending on the molecule). This however defeats the purpose of most internal representations: dimensionality reduction or faster computations.

Other (more interesting) ways of representing a molecule in internal coordinates would be:

* As (Bayati, Lardan & Bardhan 2020) focus on: a set of bond distance, angle, and **dihedrals** (see Ramachandran plots and rotamers). This is especially useful in the case of proteins. This representation is for example used by Rosetta. Angles and bond distance can be considered as constant (depending on the atomic type), and so one can infer the Cartesian coordinates of multiple atoms with just one dihedral angle, hence providing a huge dimensionality reduction. In this case there is no need for vector. Also, in the case of internal energy computation there is always one energy term depending on dihedrals, so that is already done. Limit of such methods is that in molecular constructs there are energy terms due to long-range interactions, i.e van der waals and electrostatics. There is hence a need in such software to always switch between a Dihedral representation and a Cartesian representation. Another problem is that proteins backbone are then represented as a dihedral "chain", and so the Cartesian coordinates of one atom becomes tightly linked to the previous one, which are linked to the previous one, etc. This is the problem this paper is trying to deal with.
* As (Tan et al. 2023) use, a **graph representation**. Graphein is an interesting python library for converting and analyzing proteins and molecules from cartesian coordinates to graph representation. In this case instead of having a vector, each atoms position is defined as its distance with surrounding atoms. This kind of representation has received more and more attention from the protein structure prediction & design community as it is particularly suited for neural networks. For example I would say the protein representation in AlphaFold is a kind of graph representation (using "triangles" made by the distance between three atoms). In comparison to dihedral representation, graph representation packs less "assumption" about the system and are even worst than Cartesian representation in terms of dimension (if you do a "global" representation of N atoms then you end up with N! values vs only 3\*N with Cartesian), but at least things like non-covalent interactions are implicitly represented.

I believe this answer gives enough material to look at. If this still appears unclear, skimming through Rosetta tutorials / papers figures / webinars should provide ample explanations, same for AlphaFold.

---
Tags: protein-structure
---


# Question
Title: How to remove black outline from Seurat's VlnPlot violin plot
Using Seurat's VlnPlot, how can I remove the black outline around the violin plot?

For example, how can I change from the following graph with a (black) outline:

```
    VlnPlot(ilc2,
        features = "gene",
        split.by = "group",
        cols = c("#198BFD", "#EC6464"),
        pt.size = 0) +
  geom_boxplot(width = 0.1) + theme(aspect.ratio = 1)

```

To a graph like the one below, which does not have an outline around the violin plot:

# Answer
wish there was a more subtle answer, but this works. Place this previous to your VlnPlot

```
update_geom_defaults("violin", aes(linewidth = 0))

```
> 2 votes

# Answer
```
vln_df = data.frame(gene = ilc2_bs[["RNA"]]@data[ "gene", ],
                    group = ilc2$group)

ggplot(vln_df, aes(x = group,
                   y = gene, fill = group)) +
  geom_violin(aes(fill = group,
                  color = gene),
              color = NA) +
  scale_fill_manual(values = c("#619CFF", "#F8766D")) +
  geom_boxplot(aes(fill = group),
               width = 0.1) + theme(aspect.ratio = 1)

```
> 1 votes

---
Tags: seurat, ggplot2
---