# Question
Title: Estimate insert size for single-end data with picard CollectInsertSizeMetrics
I have a BAM file generated from single-end data and I want to estimate the insert size using picard's `CollectInsertSizeMetrics` as follows:

```
java -jar Picard/picard.jar CollectInsertSizeMetrics I=Ctrl_sorted.bam O=Ctrl.insert_sizes.txt H=Ctrl.insert_sizes.pdf M=0.5 MAX_RECORDS_IN_RAM=2000000 VALIDATION_STRINGENCY=LENIENT

```

It runs fine until at the very end I get the following message:

```
WARNING 2023-01-31 15:37:13 CollectInsertSizeMetrics    All data categories were discarded because they contained < 0.5 of the total aligned paired data.

```

Picard obviously thinks the data should be paired-end, so is there a way to tell it that it is single-end? I cannot find an option for this:

```
java -jar picard.jar CollectInsertSizeMetrics -h
USAGE: CollectInsertSizeMetrics [arguments]

Collect metrics about the insert size distribution of a paired-end library. This tool provides useful metrics for
validating library construction including the insert size distribution and read orientation of paired-end
libraries.</p>The expected proportions of these metrics vary depending on the type of library preparation used,
resulting from technical differences between pair-end libraries and mate-pair libraries. For a brief primer on
paired-end sequencing and mate-pair reads, see the <a
href='https://www.broadinstitute.org/gatk/guide/article?id=6327'>GATK Dictionary</a>.<p>The CollectInsertSizeMetrics
tool outputs the percentages of read pairs in each of the three orientations (FR, RF, and TANDEM) as a histogram. In
addition, the insert size distribution is output as both a histogram (.insert_size_Histogram.pdf) and as a data table
(.insert_size_metrics.txt).</p><p>Note: Metrics labeled as percentages are actually expressed as fractions!</p><h4>Usage
example:</h4><pre>java -jar picard.jar CollectInsertSizeMetrics \<br />      I=input.bam \<br />
O=insert_size_metrics.txt \<br />      H=insert_size_histogram.pdf \<br />      M=0.5</pre>Note: If processing a small
file, set the minimum percentage option (M) to 0.5, otherwise an error may occur. <br /><br />Please see <a
href='https://broadinstitute.github.io/picard/picard-metric-definitions.html#InsertSizeMetrics'>InsertSizeMetrics</a>
for detailed explanations of each metric.<hr />
Version:2.23.8

Required Arguments:

--Histogram_FILE,-H:File      File to write insert size Histogram chart to.  Required.

--INPUT,-I:File               Input SAM or BAM file.  Required.

--OUTPUT,-O:File              The file to write the output to.  Required.

Optional Arguments:

--arguments_file:File         read one or more arguments files and add them to the command line  This argument may be
                              specified 0 or more times. Default value: null.

--ASSUME_SORTED,-AS:Boolean   If true (default), then the sort order in the header file will be ignored.  Default value:
                              true. Possible values: {true, false}

--COMPRESSION_LEVEL:Integer   Compression level for all compressed files created (e.g. BAM and VCF).  Default value: 5.

--CREATE_INDEX:Boolean        Whether to create an index when writing VCF or coordinate sorted BAM output.  Default
                              value: false. Possible values: {true, false}

--CREATE_MD5_FILE:Boolean     Whether to create an MD5 digest for any BAM or FASTQ files created.    Default value:
                              false. Possible values: {true, false}

--DEVIATIONS:Double           Generate mean, sd and plots by trimming the data down to MEDIAN +
                              DEVIATIONS*MEDIAN_ABSOLUTE_DEVIATION. This is done because insert size data typically
                              includes enough anomalous values from chimeras and other artifacts to make the mean and sd
                              grossly misleading regarding the real distribution.  Default value: 10.0.

--GA4GH_CLIENT_SECRETS:String Google Genomics API client_secrets.json file path.  Default value: client_secrets.json.

--help,-h:Boolean             display the help message  Default value: false. Possible values: {true, false}

--HISTOGRAM_WIDTH,-W:Integer  Explicitly sets the Histogram width, overriding automatic truncation of Histogram tail.
                              Also, when calculating mean and standard deviation, only bins <= Histogram_WIDTH will be
                              included.  Default value: null.

--INCLUDE_DUPLICATES:Boolean  If true, also include reads marked as duplicates in the insert size histogram.  Default
                              value: false. Possible values: {true, false}

--MAX_RECORDS_IN_RAM:Integer  When writing files that need to be sorted, this will specify the number of records stored
                              in RAM before spilling to disk. Increasing this number reduces the number of file handles
                              needed to sort the file, and increases the amount of RAM needed.  Default value: 500000.

--METRIC_ACCUMULATION_LEVEL,-LEVEL:MetricAccumulationLevel
                              The level(s) at which to accumulate metrics.    This argument may be specified 0 or more
                              times. Default value: [ALL_READS]. Possible values: {ALL_READS, SAMPLE, LIBRARY,
                              READ_GROUP}

--MIN_HISTOGRAM_WIDTH,-MW:Integer
                              Minimum width of histogram plots. In the case when the histogram would otherwise
                              betruncated to a shorter range of sizes, the MIN_HISTOGRAM_WIDTH will enforce a minimum
                              range.  Default value: null.

--MINIMUM_PCT,-M:Float        When generating the Histogram, discard any data categories (out of FR, TANDEM, RF) that
                              have fewer than this percentage of overall reads. (Range: 0 to 1).  Default value: 0.05.

--QUIET:Boolean               Whether to suppress job-summary info on System.err.  Default value: false. Possible
                              values: {true, false}

--REFERENCE_SEQUENCE,-R:File  Reference sequence file.  Default value: null.

--STOP_AFTER:Long             Stop after processing N reads, mainly for debugging.  Default value: 0.

--TMP_DIR:File                One or more directories with space available to be used by this program for temporary
                              storage of working files  This argument may be specified 0 or more times. Default value:
                              null.

--USE_JDK_DEFLATER,-use_jdk_deflater:Boolean
                              Use the JDK Deflater instead of the Intel Deflater for writing compressed output  Default
                              value: false. Possible values: {true, false}

--USE_JDK_INFLATER,-use_jdk_inflater:Boolean
                              Use the JDK Inflater instead of the Intel Inflater for reading compressed input  Default
                              value: false. Possible values: {true, false}

--VALIDATION_STRINGENCY:ValidationStringency
                              Validation stringency for all SAM files read by this program.  Setting stringency to
                              SILENT can improve performance when processing a BAM file in which variable-length data
                              (read, qualities, tags) do not otherwise need to be decoded.  Default value: STRICT.
                              Possible values: {STRICT, LENIENT, SILENT}

--VERBOSITY:LogLevel          Control verbosity of logging.  Default value: INFO. Possible values: {ERROR, WARNING,
                              INFO, DEBUG}

--version:Boolean             display the version number for this tool  Default value: false. Possible values: {true,
                              false}

Advanced Arguments:

--showHidden,-showHidden:Boolean
                              display hidden arguments  Default value: false. Possible values: {true, false}

```

Or is there another tool that can deal with single-end data to estimate the insert size?

# Answer
> 5 votes
`CollectInsertSizeMetrics` does not estimate but (as by the name) **collects** insert size metrics, which is nothing different than parsing the TLEN field from the BAM file. It is intended to only be used for paired-end data, not single-end, as single-end data do not have a value in that field in the BAM file. It is calculated by the aligner for paired-end data.

What kind of data do you have? For targeted assays such as ChIP-seq the `macs2` peak caller has a module to estimate fragment size based on strand cross correlation analysis. I do not kow whether you can use that for any sort of assay.

If this is your own data then you can ask the lab person who made the library to look up the size from the Bioanalyzer or TapeStation QC. If it is published data then check the methods text for the procedure/kit. Most NGS types have a relatively standard insert size depending on the exact assay.

---
Tags: ngs, sam, picard
---


title: Sort dictionary in descending order

text:
# Question
Title: Sort dictionary in descending order
Hi I am looking to sort a list of kmers into descending order. I used KMC to obtain a list of kmers from a fastq file separated by tab (e.g.; kmer count). From there I have written this python script to attempt to sort the kmers and their counts in descending order:

```
import csv

with open("6mers.tsv") as kmer_file:

    tsv_reader = csv.DictReader(kmer_file, delimiter="\t")

    for kmer in tsv_reader:

        name = kmer["kmer"]

        count = kmer["count"]

        # section below attempts to sort the dictionary into descending order
        # taken from (https://www.freecodecamp.org/news/sort-dictionary-by-value-in-python/)
        sorted_kmers = sorted(kmer.items(), key=lambda x:x[1], reverse=True)

        converted_dict = dict(sorted_kmers)

        print(converted_dict)

```

The code above outputs this:

```
{'kmer': 'TGCCAA', 'count': '2'}
{'kmer': 'TGGACA', 'count': '2'}
{'kmer': 'TGGGAA', 'count': '2'}
{'kmer': 'TGGTAA', 'count': '5'}
{'kmer': 'TGTCAA', 'count': '5'}
{'kmer': 'TGTGAA', 'count': '2'}
{'kmer': 'TTCCAA', 'count': '13'}
{'kmer': 'TTGAAA', 'count': '6'}
{'kmer': 'TTGCAA', 'count': '4'}
{'kmer': 'TTTAAA', 'count': '6'}

```

Clearly I am creating a separate dictionary for each kmer and it's count - I'm not sure how but I would like a result that looks more like this:

```
{'TTCCAA', '13'}
{'TTGAAA', '6'}
{'TTTAAA', '6'}
{'TGGTAA', '5'}
{'TGTCAA', '5'}
{'TTGCAA', '4'}
{'TGCCAA', '2'}
{'TGGACA', '2'}
{'TGGGAA', '2'}
{'TGTGAA', '2'}

```

Can someone help me fix my code so it works?

**EDIT:**

The .tsv file looks just like this:

```
TGCCAA  2
TGGACA  2
TGGGAA  2
TGGTAA  5
TGTCAA  5
TGTGAA  2
TTCCAA  13
TTGAAA  6
TTGCAA  4
TTTAAA  6

```

# Answer
> 2 votes
Just make sure to cast the count to an integer prior to sorting. Assuming you need to do something with the counts, you could use:

```
import csv

from operator import itemgetter
from pprint import pprint

with open("6mers.tsv") as kmer_file:

    reader = csv.reader(kmer_file, delimiter="\t")
    counts = ((kmer, int(count)) for kmer, count in reader)

    sorted_kmers = sorted(counts, key=itemgetter(1), reverse=True)

    pprint(sorted_kmers)

```

Results:

```
[('TTCCAA', 13),
 ('TTGAAA', 6),
 ('TTTAAA', 6),
 ('TGGTAA', 5),
 ('TGTCAA', 5),
 ('TTGCAA', 4),
 ('TGCCAA', 2),
 ('TGGACA', 2),
 ('TGGGAA', 2),
 ('TGTGAA', 2)]

```

# Answer
> 1 votes
This is untested code ... without seeing to the raw file its tricky. Could you try:

```
import csv
with open("6mers.tsv") as kmer_file:
    tsv_reader = csv.DictReader(kmer_file, delimiter="\t")
    mydict = {rows[0]:rows[1] for rows in tsv_reader if not rows[0] == 'kmer'}

sorted_kmers = sorted(mydict.items(), key=lambda x:x[1], reverse=True)
print(sorted_kmers)

```

I'm assuming there is a header called 'kmer'

---

The output didn't work. The rows\['kmers'\] infers that (somehow) tsv\_reader is already a dictionary (that's not how I thought it worked). Thus try,

```
import csv
with open("6mers.tsv") as kmer_file:
    tsv_reader = csv.DictReader(kmer_file, delimiter="\t")
    mydict = {rows['kmers']:rows['count'] for rows in tsv_reader}

sorted_kmers = sorted(mydict.items(), key=lambda x:x[1], reverse=True)
print(sorted_kmers)

```

Unless it's loaded as a singular dictionary with Kmer as the key and count as the value the sort line will not work. Once that format is achieved as a dictionary then it will work.

What I don't understand is where it's getting 'kmer' or 'count' from if there are no headers in the csv file. Something isn't right if the above code works.

---
Tags: python, k-mer, dictionary, tsv
---


# Question
Title: How do I cut a genome into one specific region?
Good morning, I study bioinformatics and I need some help with a work. I have the NCBI accession numbers for the species I want for example: NC\_037078.1 and I have the specific coordinates where I need to cut them to obtain the 16s gene but I don´t know how to do it. I'm using the language shell script. Is there some way to cut the genome into one specific region?

# Answer
If you know the coordinates, you could just use `samtools faidx` to extract the corresponding subsequence from the FASTA file(s). Regions can be specified on the command line with the format: `RNAME[:STARTPOS[-ENDPOS]]`. For example, to extract the first ten bases from `chr1` in your reference FASTA file, you would use:

```
samtools faidx ref.fasta chr1:1-10

```

Or if you wanted the whole chromosome, you would use:

```
samtools faidx ref.fasta chr1

```

---

---

With regards to your specific example, you could try:

```
samtools faidx NC_037078.1.fasta NC_037078.1:96630-98120

```

Results:

```
>NC_037078.1:96630-98120
TCTCATGGAGAGTTCGATCCTGGCTCAGGATGAACGCTGGCGGCATGCTTAACACATGCA
AGTCGGACGGGAAGTGGTGTTTCCAGTGGCGGACGGGTGAGTAACGCGTAAGAACCTGCC
CTTGGGAGGGGAACAACAGCTGGAAACGGCTGCTAATACCCCGTAGGCTGAGGAGCAAAA
GGAGGAATCCGCCCGAGGAGGGGCTTGCGTCTGATTAGCTAGTTGGTGAGGCAATAGCTT
ACCAAGGCGATGATCAGTAGCTGGTCCGAGAGGATGATCAGCCACACTGGGACTGAGACA
CGGCCCAGACTCCTACGGGAGGCAGCAGTGGGGAATTTTCCGCAATGGGCGAAAGCCTGA
CGGAGCAATGCCGCGTGGAGGTAGAAGGCCCACGGGTCGTGAACTTCTTTTCCCGGAGAA
GAAGCAATGACGGTATCTGGGGAATAAGCATCGGCTAACTCTGTGCCAGCAGCCGCGGTA
AGACAGAGGATGCAAGCGTTATCCGGAATGATTGGGCGTAAAGCGTCTGTAGGTGGCTTT
TTAAGTCCGCCGTCAAATCCCAGGGCTCAACCCTGGACAGGCGGTGGAAACTACCAAGCT
GGAGTACGGTAGGGGCAGAGGGAATTTCCGGTGGAGCGGTGAAATGCGTAGAGATCGGAA
AGAACACCAACGGCGAAAGCACTCTGCTGGGCCGACACTGACACTGAGAGACGAAAGCTA
GGGGAGCGAATGGGATTAGATACCCCAGTAGTCCTAGCCGTAAACGATGGATACTGGGCG
CTGTGCGTATCGACCCGTGCAGTGCTGTAGCTAACGCGTTAAGTATCCCGCCTGGGGAGT
ACGTTCGCAAGAATGAAACTCAAAGGAATTGACGGGGGCCCGCACAAGCGGTGGAGCATG
TGGTTTAATTCGATGCAAAGCGAAGAACCTTACCAGGGCTTGACATGCCGCGAATCCTTT
TGAAAGAGAGGGGTGCCTTCGGGAACGCGGACACAGGTGGTGCATGGCTGTCGTCAGCTC
GTGCCGTAAGGTGTTGGGTTAAGTCCCGCAACGAGCGCAACCCTCGTGTTTAGTTGCCAT
CGTTGAATTTGGAACCCTGAACAGACTGCCGGTGATAAGCCGGAGGAAGGTGAGGATGAC
GTCAAGTCATCATGCCCCTTATGCCCTGGGCGACACACGTGCTACAATGGCCGGGACAAA
GGGTCGCGATCTCGCGAGGGTGAGCTAACCCCAAAAACCCGTCCTCAGTTCGGATTGCAG
GCTGCAACTCGCCTGCATGAAGCCGGAATCGCTAGTAATCGCCGGTCAGCCATACGGCGG
TGAATTCGTTCCCGGGCCTTGTACACACCGCCCGTCACACTATGGGAGCTGGCCATGCCC
GAAGTCGTTACCTTAACCGCAAGGAGGGGGATGCCGAAGGCAGGGCTAGTGACTGGAGTG
AAGTCGTAACAAGGTAGCCGTACTGGAAGGTGCGGCTGGATCACCTCCTTT

```

Note that one or more regions may be specified on the command line and all position coordinates are 1-based. If the index file (i.e. `<ref.fasta>.fai`) doesn't exist, it will be created. If no regions are specified, `samtools faidx` will simply create the index file and exit. This is useful if you need to create the index prior to querying it. Only uncompressed FASTA files or FASTA files compressed with bgzip can be indexed.

```
$ samtools faidx -h
Usage: samtools faidx <file.fa|file.fa.gz> [<reg> [...]]
Option:
 -o, --output FILE        Write FASTA to file.
 -n, --length INT         Length of FASTA sequence line. [60]
 -c, --continue           Continue after trying to retrieve missing region.
 -r, --region-file FILE   File of regions.  Format is chr:from-to. One per line.
 -i, --reverse-complement Reverse complement sequences.
     --mark-strand TYPE   Add strand indicator to sequence name
                          TYPE = rc   for /rc on negative strand (default)
                                 no   for no strand indicator
                                 sign for (+) / (-)
                                 custom,<pos>,<neg> for custom indicator
     --fai-idx      FILE  name of the index file (default file.fa.fai).
     --gzi-idx      FILE  name of compressed file index (default file.fa.gz.gzi).
 -f, --fastq              File and index in FASTQ format.
 -h, --help               This message.

```
> 3 votes

---
Tags: shell
---


title: How can I tranform this String?

text:
# Question
Title: How can I tranform this String?
I'm trying to alter this string: `NC_037078 1:122545 124035` into this: `NC_037078.1:122545-124035` but I don´t know how.

The document is like this:\

```
    1. rrn16 \
16S ribosomal RNA [Genlisea aurea]\
Other Aliases: C7J69_pgr001, BK572_Gau_0120\
Genomic context: Chloroplast\
Annotation:  NC_037078.1 (122545..124035, complement)\
ID:` 36166080

```

The script until now looks like this:

```
code=$(ls -l | grep "Annotation" data.txt | cut -d " " -f3,4 | sed 's/[(]/:/g' | sed 's/[)]/ /g' | sed 's/[,]/ /g' | sed 's/[ ]//g' | sed 's/["."]/ /g' )

```

Can you help me solving this problem?

# Answer
`echo NC_037078 1:122545 124035 | sed 's/ /./1' | sed 's/ /:/g'`

Note: Reference: https://stackoverflow.com/questions/59047776/sed-replace-only-first-occurrence-then-only-second-occurrence-then-only-third
> 1 votes

---
Tags: shell
---


title: How can I transform this String?

text:
# Question
Title: How can I transform this String?
Good afternoon, I'm a bioinformatics student and I'm trying to put some strings between commas with a space. The output looks like this: `NC_037078.1:122545-124035,NC_037078.1:96630-98120` and I want it to be like this: `NC_037078.1:122545-124035, NC_037078.1:96630-98120`

The document is like this:
1. rrn16
16S ribosomal RNA \[Genlisea aurea\]
Other Aliases: C7J69\_pgr001, BK572\_Gau\_0120
Genomic context: Chloroplast
Annotation: NC\_037078.1 (122545..124035, complement)
ID:\` 36166080

The script until now looks like this:

```
code=$(ls -l | grep "Annotation" data.txt | cut -d " " -f3,4 | sed 's/[(]/:/g' | sed 's/[)]/ /g' | sed 's/[,]/ /g' | sed 's/[ ]//g' | sed 's/["."]/ /g' | sed 's/ /./1' | sed 's/ /-/g' | sed 's/-//1' | tr "\n" ",")

```

# Answer
It seems a small change. Anyway try ...

```
perl -p -e 's/(.*NC_[0-9]+\.[0-9]),(NC_[0-9]+\.[0-9].*)/$1,\s$2/g' data.txt

```

If you have happy that this works, then repeat with an `-i` flag

```
perl -p -i -e 's/(.*NC_[0-9]+\.[0-9]),(NC_[0-9]+\.[0-9].*)/$1,\s$2/g' data.txt

```

I've not checked the code, but it's worth checking whether `.*` is really needed. Its a while since I've used Perl
> 1 votes

---
Tags: shell
---


# Question
Title: hierarchical clustering of kmers and their counts
I have a list of kmers and their frequencies broken down like this:

```
TGCCAA  2
TGGACA  2
TGGGAA  2
TGGTAA  5
TGTCAA  5
TGTGAA  2
TTCCAA  13
TTGAAA  6
TTGCAA  4
TTTAAA  6

```

And I would like to use clustering to look at the most frequent and the most unique instances of kmers in a sample. I started off following the tutorial found here (https://www.w3schools.com/python/python\_ml\_hierarchial\_clustering.asp#:~:text=Hierarchical%20clustering%20is%20an%20unsupervised,need%20a%20%22target%22%20variable.) using the scikit-learn python package.

However, when I go to implement my data I get a ValueError where it can not convert string to float, which seems obvious as the tutorial uses two sets of numbers and not strings.

I was wondering if there was a known way to cluster kmers and their counts or simply just strings.

# Answer
> 2 votes
This isn't good data for clustering based on the count value because it needs to be multivariate, e.g. you'd need two or more 'count' columns (if that makes any biological sense). This data is univariate so I'd just use a standard bar chart, personally I would rank the data and consider the underlying distribution.

You could make a tree of the sequence data and the data counts are the labels. If you want to make a tree (hierarchical clustering) then using nucleotide data its becomes a phylogenetic tree, which isn't hierarchical clustering - although clustering can be used to make phylogenetic trees.

---

K-mers ... if these sites are not homologous a phylogenetic tree is not a good solution: they do look like homologous however. In terms of shifting to numerical variables the technical term is called `one-hot-key` encoding and scikit-learn has a function for this. This is what Steve is saying. However, understanding what the k-mer represents is kinda important. To perform one-hot-key encoding its:

```
from sklearn.preprocessing import OneHotEncoder
hot = OneHotEncoder()
.... gets complicated

```

Alternative approach

```
import pandas as pd
df = pd.read_csv("kmer.csv")
one_hot_encode = pd.get_dummies(df, columns = ['kmers'])

```

This code will transform your data to what Steve said .. keeping it as `pandas` is sort of useful because hierarchical clustering will take a pandas input.

I do have doubts ... the `counts` column needs a unique identifier *in addition* to the number, because you'll get identical labels otherwise (my original idea had limitations). Ultimately, I don't precisely know what the kmer value biologically means and thats essential. I would be cautious about dumping that straight into a hierarchical clustering algorithm (it will be in `scipy` BTW not `scikit-learn`) .

---

The following would provide a phylogenetic tree input, but statistically its dodgy without understanding what the k-mer means

```
import re

kmer = ['TGCCAA', 'TGGACA', 'TGGGAA', 'TGGTAA', 'TGTCAA', 'TGTGAA',  'TTCCAA', 'TTGAAA', 'TTGCAA', 'TTTAAA']

count = [2, 2,  2, 5, 5, 2, 13,  6,  4,  6]

x = [chr(i) for i in range(ord('a'),ord('z')+1)]
n = 0
mydict = {}
for c,k in zip (count, kmer):
    n += 1
    mydict[str(c)+x[n]] = k
    if n > 25:
        n = 0
for key in mydict:
    if mydict[key] == None:
        continue
    print ('>' + key + '\n' + mydict[key] + '\n')

```

output

```
>2b
TGCCAA
>2c
TGGACA
>2d
TGGGAA
>5e
TGGTAA
>5f
TGTCAA
>2g
TGTGAA
>13h
TTCCAA
>6i
TTGAAA
>4j
TTGCAA
>6k
TTTAAA

```

Here's the tree

All I can conclude is count and kmer - with one exception `TTTAAA` and `TTGAAA` are not correlated. `TTRAAA` maybe be a result, possibly.

Here's the matrix.

---

**EDIT** Actually there might be a result here: the k-mer representing '13' is a genetic outlier. Whilst it might be difficult to see in the tree I've presented that could proven more decisively. If a bigger sample size demonstrated comparable phenomena - that is a result for certain. So maybe it's not such a bad analysis after all - it could lead to notably interesting result. Hierarchical clustering BTW couldn't generate that results - its assumption base removes the outlier effect.

# Answer
> 3 votes
I think you need some way to convert the k-mer strings into a numerical representation so that they can be clustered. There's a few ways to do this. For example, you could one-hot (or two-bit) encode each k-mer so that each nucleotide is represented as a set of binary values. Or you could just use a simple integer encoding. For example:

```
def int_encode(kmer):
    two_bit = {'A': 0, 'C': 1, 'G': 2, 'T': 3}
    encoding = []
    for nucleotide in kmer:
        encoding.append(two_bit[nucleotide])
    return encoding

```

You could then concatenate the output with the frequency counts of each k-mer. For your set of k-mers, your input for clustering would then look like:

```
[3, 2, 1, 1, 0, 0, 2]
[3, 2, 2, 0, 1, 0, 2]
[3, 2, 2, 2, 0, 0, 2]
[3, 2, 2, 3, 0, 0, 5]
[3, 2, 3, 1, 0, 0, 5]
[3, 2, 3, 2, 0, 0, 2]
[3, 3, 1, 1, 0, 0, 13]
[3, 3, 2, 0, 0, 0, 6]
[3, 3, 2, 1, 0, 0, 4]
[3, 3, 3, 0, 0, 0, 6]

```

---
Tags: python, clustering, k-mer
---


# Question
Title: Heat map of protein expression from normalized abundance
This is a follow up to my question posted here Processing proteomics data

In relevance to my research, I've been looking for proteomics data (control vs. diabetic) and I found a dataset in the article "Diabetes causes marked inhibition of mitochondrial metabolism in pancreatic β-cells" (figures from here: link).

I would like to know how to generate the following heat map from the normalized abundances reported.

Fold change heatmap:

Abundance:

Figure description:

> Heat maps of relative mRNA and protein levels of the indicated genes in islets isolated from control and 2-week diabetic βV59M mice. Each box corresponds to a different animal. Colour indicates log2 fold-change. Heat map: TCA cycle. Red, proteins upregulated in diabetes; blue, proteins downregulated in diabetes. Grey, no change or not detected. PDK1, pyruvate dehydrogenase kinase 1. Abundance: Abundance of the indicated proteins, quantified by mass spectrometry, in islets isolated from control (black, Ctrl, n = 4) and 2-week diabetic βV59M (white, Diab, n = 4) mice.

I am not able to understand how the log2 fold-change is computed from the abundance. For instance, if we consider `PDK1` (first panel) the abundance measure of control is ~50000 and diabetic is ~200000. Heat map says log2 fold-change (in the figure description added above). What is the fold-change that is referred to? From what is shown, I think this is not diabetic/control fold change; we see two panels (i.e one for control (4 samples) and the other for diabetic (4 samples)) and not a single panel. Since the legend is from (-2 to 2) I think these values are not log2 (abundance) as well.

Could someone please clarify how the heat map has been generated from the abundance measures? Basically, I want to know how to generate the first heatmap figure from abundance plots.

# Answer
The -2 to 2 should be the normalized version of log2 data (maybe mean or median normalization). This is commonly used in heatmaps to show downregulated value vs. upregulated value.
> 1 votes

---
Tags: proteins, proteomics, heatmap, mass-spectrometry, meta-analysis
---


# Question
Title: what does ALT ID DUP refer to when ALT ID for copy number are already specified?
ALT ID CN2, CN3,CN4,CN5 are regions of elevated copy number relative to the reference, Then what does the ID =DUP separately cover? I am new to these analyses. please help

```
##ALT=<ID=DUP, Description="Region of elevated copy number relative to the reference">
##ALT=<ID=CNO, Description="Copy number allele: 0 copies">
##ALT=<ID=CN2, Description="Copy number allele: 2 copies">
##ALT=<ID=CN3, Description="Copy number allele: 3 copies">
##ALT=<ID=CN4, Description="Copy number allele: 4 copies">
##ALT=<ID=CN5, Description="Copy number allele: 5 copies">

```

# Answer
> 1 votes
The details will depend on the specific tool used to generate the VCF. Each tool has its own system and there is not much standardization for these fields. That said, the likeliest situation here is that if the CNV caller cannot determine the copy number, but only that there is a duplication, it just puts DUP. And if it can determine the specific copy number, then it uses one of the specific terms.

Alternatively, if the copy number is beyond those numbers with a specific `CNN` value (say you have 7 copies, or 65 or whatever), then it would presumably also use a `DUP`.

---
Tags: vcf, illumina, cnv
---


# Question
Title: relation between Illumina sequencing primer and viral sequences
Dealing whith a problematic sequencing run I found this over-represented sequence:

> **GGAAGAGCACACGTCTGAACTCCAGTCAC**TAGCTT**ATCTCGTATGGCGTCTTCTGCTTG**

It is clearly relate to the Illumina sequencing primer as shown by this read structure description (note the bold residues above):

I searched the sequence in the nucleotide NR ncbi database, using blastn (standard parameters) and I found this as a first hit:

```
Severe acute respiratory syndrome coronavirus 2 genome assembly, chromosome: 1
Sequence ID: OV196492.1Length: 30243Number of Matches: 2
Range 1: 73 to 133GenBankGraphicsNext MatchPrevious Match
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
89.8 bits(48)   2e-14   57/61(93%)  2/61(3%)    Plus/Minus
Query  1    GGAAGAGCACACGTCTGAACTCCAGTCACTAG--CTTATCTCGTATGGCGTCTTCTGCTT  58
            ||||||||||||||||||||||||||||| ||  ||||||||||||| ||||||||||||
Sbjct  133  GGAAGAGCACACGTCTGAACTCCAGTCACCAGTGCTTATCTCGTATGCCGTCTTCTGCTT  74

Query  59   G  59
            |
Sbjct  73   G  73

Range 2: 188 to 248GenBankGraphicsNext MatchPrevious MatchFirst Match
Alignment statistics for match #2
Score   Expect  Identities  Gaps    Strand
89.8 bits(48)   2e-14   57/61(93%)  2/61(3%)    Plus/Minus
Query  1    GGAAGAGCACACGTCTGAACTCCAGTCACTAG--CTTATCTCGTATGGCGTCTTCTGCTT  58
            ||||||||||||||||||||||||||||| ||  ||||||||||||| ||||||||||||
Sbjct  248  GGAAGAGCACACGTCTGAACTCCAGTCACCAGTGCTTATCTCGTATGCCGTCTTCTGCTT  189

Query  59   G  59
            |
Sbjct  188  G  188

```

There are also sequences of other virues that match very well.

Do you know what is the relation between Illumina adapters/primers and viral seqiuences?

# Answer
> 8 votes
You'll unfortunately find adapter sequences contaminating a lot of entries in Genbank. In short, the sequence is most likely either an unassembled fragment where the input into the assembly step hadn't had adapters trimmed or (somewhat more rarely) the adapter sequence was incorporated into an assembled contig and then uploaded. You'll see a lot of lower quality sequences like this, especially from massive viral screens.

Be aware that there are other quality issues facing particularly virus sequences such as extremely high NNN stretches (cf. the difficulties faced in creating the C-RVDB viral database once there were hundreds of thousands of apparently poor-quality coronavirus sequences automatically added to Genbank).

# Answer
> 0 votes
SARS-CoV-2 NCBI no. OV196492.1

**Summary** The relationship sought is simply the genome position because this is the UTR (UnTranslated Region).

The linker is present starting at position 74 and is also present at position of 189 of the SARS-CoV-2 genome from 5' to 3' (I checked). This is where position 1 is not necessarily position 1 of the genome, but the start of the authors genome. So it's the UTR. I did suspect this was the location.

There are three reasons for the drop in sequence quality at the UTR:

1. overall quality, really the read depth, of short-read sequence data is reduced at the terminal regions of viruses *per se*;
2. there is extensive secondary structure in the UTR regions which is critical to its function, thus the folding will affect short-read sequencing approaches, e.g. primer binding will be more difficult because these sequences will preferentially self-anneal (preventing the primer from binding);
3. QC annotation from the authors will focus on the coding regions, the UTR is outside the coding regions and therefore can get missed.

In this particular case point 3 is likely the underlying reason. The authors QC ignored the UTRs (the 3' UTR is also problematic). If a linker is caught in the coding domains the translation will flag it immediately ... obviously the UTR isn't translated. Thats simply why. For virologists, this isn't cool, but for clinicians they don't see its importance.

**Overall**, this is simply an issue of *ad hoc* QC because the NCBI SARS-CoV-2 database is extremely powerful and excellent for custom analysis. To "throw in the towel" on the power and flexibility of NCBI's SARS-CoV-2 database for what is essentially some personal QC coding isn't worth the trade-off in my personal opinion. It's a powerful portal, which enables an analysts to focus on downstream issues.

---
Tags: covid-19, illumina, virus, coronavirus, adapter
---


# Question
Title: How to get a file with the number of reads for several fastq.gz files?
I have generated several FASTQ files and I would like to know the amount of reads for each of them. I am planning to run FastQC on the files which I know would give me the number of reads per sample but it would be a lot easier for me if I could get one text file with the file name and number of reads for all samples together.

When I check for one file I normally do `zcat filename_R1.fastq.gz | wc -l` and then divide by 4. Is there a way to do it for all my files at once and get a text file with the file name and reads at the end?

For example I will have 71 R1.fastq.gz files and I would like to output the amount of reads to a text file. I just meant that I would prefer to use one command instead of having to use the command I mentioned above for each file.

# Answer
One way would be to use `seqkit stats`, which can be installed using `conda`:

```
conda install -c bioconda seqkit

```

It has a number of useful options, including ways to parallelize the counting, using `-j`, `--threads` (default 4), and to skip errors using `-e`, `--skip-err`. The latter can be especially useful if you have a larger number of FASTQ files with millions of reads and you'd like to run the command overnight for example. There's also an option to produce a machine-friendly tabular output, using `-T`, `--tabular` which would allow you to parse the output without much effort. I'd be tempted to just write the output to a file using the `-o`, `--out-file` option:

```
seqkit stats -To stats.tsv *.fastq.gz

```

Example results:

```
$ cat stats.tsv
file    format  type    num_seqs    sum_len min_len avg_len max_len
PHB_Rep1_R1.fastq.gz    FASTQ   RNA 118571  11857100    100 100.0   100
PHB_Rep1_R2.fastq.gz    FASTQ   RNA 118571  11857100    100 100.0   100
PHB_Rep2_R1.fastq.gz    FASTQ   RNA 144826  14482600    100 100.0   100
PHB_Rep2_R2.fastq.gz    FASTQ   RNA 144826  14482600    100 100.0   100
PHB_Rep3_R1.fastq.gz    FASTQ   RNA 129786  12978600    100 100.0   100
PHB_Rep3_R2.fastq.gz    FASTQ   RNA 129786  12978600    100 100.0   100
PUH_Rep1_R1.fastq.gz    FASTQ   RNA 227392  22739200    100 100.0   100
PUH_Rep1_R2.fastq.gz    FASTQ   RNA 227392  22739200    100 100.0   100
PUH_Rep2_R1.fastq.gz    FASTQ   RNA 162373  16237300    100 100.0   100
PUH_Rep2_R2.fastq.gz    FASTQ   RNA 162373  16237300    100 100.0   100
PUH_Rep3_R1.fastq.gz    FASTQ   RNA 185442  18544200    100 100.0   100
PUH_Rep3_R2.fastq.gz    FASTQ   RNA 185442  18544200    100 100.0   100

```

And then just parse the output to get what you need, if or when you need it. For your specific requirements, you could just use AWK to skip the header line and print the first and fourth columns:

```
awk 'BEGIN {FS=OFS="\t" } NR>1 { print $1, $4 }' stats.tsv

```

Results:

```
PHB_Rep1_R1.fastq.gz    118571
PHB_Rep1_R2.fastq.gz    118571
PHB_Rep2_R1.fastq.gz    144826
PHB_Rep2_R2.fastq.gz    144826
PHB_Rep3_R1.fastq.gz    129786
PHB_Rep3_R2.fastq.gz    129786
PUH_Rep1_R1.fastq.gz    227392
PUH_Rep1_R2.fastq.gz    227392
PUH_Rep2_R1.fastq.gz    162373
PUH_Rep2_R2.fastq.gz    162373
PUH_Rep3_R1.fastq.gz    185442
PUH_Rep3_R2.fastq.gz    185442

```
> 5 votes

# Answer
You can replace `gunzip -c` with `zcat`.

```
#USAGE: script.sh /runFolder
for i in $1/*.1.fastq.gz
do
   echo "$i" $(gunzip -c $i | echo `wc -l`/4 | bc -l)
done

#OR in the current folder, do this:
for i in ./*.1.fastq.gz
do
   echo "$i" $(gunzip -c $i | echo `wc -l`/4 | bc -l)
done

```
> 1 votes

# Answer
`zcat`, like `cat`, works with multiple file names as arguments, up to the argument limit (on Linux, that's usually at least a thousand), so counting multiple lines in multiple files still works, assuming the fastq files are formatted in the expected way with only four lines per record:

```
zcat filename_R1.fastq.gz filename_R2.fastq.gz ... | wc -l

```

If there's a consistent file pattern, that can be used with file name glob wildcards:

```
zcat *.fastq.gz | wc -l

```

If you have a file containing a list of files, this can be processed as shell arguments as well:

```
zcat $(cat list_of_files.txt) | wc -l

```

For an unlimited number of files, something like GNU parallel can be used to feed file names from piped input. It's also a bit safer for file names containing weird characters in them:

```
cat list_of_files.txt | parallel zcat {} | wc -l

```

... but \[after checking the question again\] given that you want counts for individual files, the file name needs to be added as well:

```
cat list_of_files.txt | parallel echo {} $(zcat {} | wc -l)

```

The `read` function provides an alternate way to do the same thing:

```
cat list_of_files.txt | while read fileName; do echo ${fileName} $(zcat ${fileName} | wc -l); done

```

Assuming you are using bash, dividing by 4 can be done via the arithmetic substitution operator, `$((X))`:

```
cat list_of_files.txt | parallel echo {} $(( $(zcat {} | wc -l) / 4 ))

```

Or with `read`:

```
cat list_of_files.txt | while read fileName; do echo ${fileName} $(( $(zcat ${fileName} | wc -l) / 4 )); done

```

And, finally, this can be output to a file using `>`:

```
(cat list_of_files.txt | while read fileName; do echo ${fileName} $(( $(zcat ${fileName} | wc -l) / 4 )); done) > fastq_record_counts.txt

```
> 0 votes

# Answer
If you're looking for a one-liner,

```
for i in $(ls *.fastq.gz); do count=$(zcat $i | grep "^@" | wc -l); printf "%b\t%b" $i $count; done > counts.tsv

```
> -4 votes

---
Tags: fastq, linux, command-line
---


# Question
Title: How to do pathway analysis after scanpy for single cell data after DE analysis?
I have done DE analysis using SCANpy on my single cell data and I have compared each cluster versus all the other clusters. One cluster seems particularly interesting so I wanted to do pathway analysis.

I have already tried to do the pathway analysis but I got confused. I initially followed the instructions from here (https://nbisweden.github.io/workshop-scRNAseq/labs/compiled/scanpy/scanpy\_05\_dge.html) and did the enrichment and then I also wanted to do GSEA so I used gseapy.prerank. However, then I was trying to find a way to visualize the network so I ended up going to the gseapy documentation.

The steps of the enrichment analysis seemed the same but then I got confused with the GSEA portion. The first tutorial I found used gseapy.prerank but the second one uses gp.gsea. I was trying to understand the difference and it seems that since I already have the DE analysis I should be using the gseapy.prerank right bit I am not sure.

I was wondering if doing the analysis trough KEGG would be easier.

Thank you

> Gene set enrichment analysis (GSEA) is a commonly used algorithm for characterizing gene expression changes. However, the currently available tools used to perform GSEA have a limited ability to analyze large datasets, which is particularly problematic for the analysis of single-cell data. To overcome this limitation, we developed a GSEA package in Python (GSEApy), which could efficiently analyze large single-cell datasets.

# Answer
> 1 votes
An easy way to do an enrichment analysis in scanpy is using `scanpy.queries.enrich()`. This relies on gprofiler, which uses KEGG, GO and other ontologies

# Answer
> 0 votes
The original question was GSEApy vs KEGG. The response below is GSEApy. The question has been edited and is now a different question concerning options within GSEApy.

---

*Response to original question*

If GSEApy as described in the publication provides the output you require - that is much, much simpler. KEGG will be complex and you'd need its ids.

However, surely GSEA will require to perform the full calculation of fold-change from scratch rather than be fed the output of scanpy? KEGG will require the output of scanpy. GSEApy primary function is fold-change rather identifying the precisely which pathways are regulating the phenotype under analysis.

KEGG will provide a more accurate answer - in my humble opinion. GSEApy is a lot easier because it very much an "all-in-one" solution. If the answer isn't satisfactory you can always switch.

---

Just to note the GSEApy version I'm referring to is the latest publication.

**NOTE** There was a deleted post from the OP about an option within GSEA. This is a separate question in the SE format of questions and answers. Thats the way the site works.

---
Tags: single-cell, pathway, gsea
---


# Question
Title: How do cellranger and cellbender call cells? What is the difference between them?
When I use cellranger or cellbender to filter cell from raw RNA count matrix, the results they output are allways different. Each of then can call some cells not appearing in other's result. SO, how do cellranger and cellbender call cells? What is the difference between them?

# Answer
As mentioned in the comments, an example of your inputs, outputs and their difference would be a good addition to your question. Moreover, the cell calling changed in CellRanger v3 compared to V2.

There are several differences in how CellRanger and CellBender filter out empty droplets and call cells. Both methods are either published or available as preprints, so you can learn much more about them especially if you're interested in the technical aspects.

CellRanger v3 uses a procedure largely based on the emptyDrops algorithm (Lun et al. 2019).

Briefly, the algorithm calculates an "ambient profile" by summing all the read counts belonging to barcodes with a total read count below a threshold. Then, it compares each barcode to this ambient profile, testing for the similarity of the distribution of a barcode's counts to the ambient's. CellRanger v3 does not implement exactly the same thresholds and will give you slightly different results than emptyDrops. CellRanger v2, on the other hand, used hard thresholds (based on knee point or the inflection point of the barcode ranks) which is less precise and biologically justifiable.

CellBender (Fleming et al. 2022) uses a different procedure, in that it uses a deep learning model to learn suitable priors for ambient count distribution, and uses a different statistical model to call cells. The authors of CellBender compare their method to CellRanger and emptyDrops, so you can see for yourself how different the results are (on control datasets).

Additionally, one of the differences is that CellBender will give you "corrected counts", i.e. it will subtract the estimated ambient profile from the "real" profile. In other words, it will remove the background noise and give you counts that are (supposedly) only coming from good cells.

Just keep in mind that whether or not "corrected counts" are good to use for downstream applications (HVG identification, clustering, etc.) is debatable (see here for a discussion).

A reasonable compromise is to flag genes with high ambient contaminations so that you know whether to trust the result of e.g. a differential expression test or a marker finding procedure.
> 1 votes

---
Tags: scrnaseq, cellranger
---


# Question
Title: Cells with zero expression of a given gene
I have scRNA-seq from PBMC analyzed by Seurat. I am seeing for most of genes a lot of cells have a zero expression like this gene

So 268 single cells have the expression level == 0

And I am more wondering when seeing zero expression for a housekeeping gene like

For which 72 single cells do not express this gene

Does this come from biology or technical artifact?

Thanks for sharing your ideas here

# Answer
> 6 votes
This has been debated for a few years now as the "dropout" problem, which is actually a mixture of different issues.

On one hand, you are sequencing a relatively low input at a depth that does not reach saturation (i.e. you don't measure all there is to measure). This is partially because, on droplet data - such as 10X - you never use the same sequencing depth per single cell that you would use for a bulk measurement. You produce fewer sequencing reads.

Additionally, the way a library is generated within a droplet with such a low input implies that transcripts will compete for a relatively limited amount of enzymes and reagents.

This means that a highly expressed gene (i.e. one of which a single cell has produced several transcript molecules) will be more likely to be detected than one with fewer molecules (there are some biases but they do not explain alone the phenomenon). This would hardly ever happen in a bulk RNA sequencing dataset produced at an appropriate depth, and it is why we regard count data from scRNA-seq as "sparse" (with many 0's).

On the other hand, there is a third issue to consider that is more "biological" than 'technical": the way we used to think about transcription was rather mechanistic (here I hope others will excuse this word).

What I mean by that is that you would think of transcription as an ordered set of steps whereby a transcription factor binds a promoter aided by other cofactors and chromatin complexes, recruit the RNA polymerase and let it do its job for as long/as many molecules as needed. However more biophysically motivated models generated from both *in vitro* and *in vivo* lines of evidence have painted a different picture, in which transcription is bursty, has restarts, early stops, and does not necessarily work as a putting a needle on a vinyl record. DNA loci may be transcribed at low levels due to phenomena that don't immediately depend on transcriptional regulation, such as cell cycle dynamics. This is a big source of intercellular variability. If you think about it, these are large molecules operating in solution and associating/dissociating from other molecules according to statistical thermodynamics. If the "vinyl" view of transcription was good enough to make sense of most bulk RNA-seq results, where you are indeed observing the average expression of genes in a cell population/tissue, once we get at the single cell level such a model does not necessarily hold.

You know that GAPDH is a "housekeeping" gene because on a population average measured by RT-qPCR or Western Blot you never\* see its levels fluctuate significantly, and they stay high compared to many other genes. At the single cell level, however, this need not be true at all times and for all cells. Hence a proportion of your zero-count genes will come from detection limits and low input ("technical" dropouts) and another proportion will come from the stochastic nature of transcription itself.

Now the important aspect that the field has reflected upon is how to model these dropout, or if we even need to. The goal of these models is to understand variability and noise in the signal, allowing us to extrapolate useful biological data. Should we use a specific distribution of counts that accounts for an excess of 0's in the data (zero inflation)? Should we model dropouts entirely separately? Are distributions used in bulk RNA-seq analysis sufficient to model the data in an actionable and computationally efficient way? What about modelling burstiness and other biophysical properties?

I believe it can be argued that there is no "actual/correct statistical model" because models are choices and results should be interpreted in light of those choices and their assumptions. If you are looking to cluster your cells and find "markers" to characterize the composition of a tissue, a simple model may be sufficiently equipped to let you draw some interesting conclusions. If instead you want to look more in depth at dynamical properties (transitions between cell states, "RNA velocity", splicing, etc.) you may need to use more rigorous, biophysically motivated models that are also more complex in their theoretical justifications and their implementation.

Here are some recommended readings if you want to dive deeper in the issue of dropouts and count modelling. The list is not exhaustive but I hope it helps.

\*more precise measurements in different tissues/cell lines have shown housekeeping genes to be variable in bulk as well, but that's another issue.

---
Tags: scrnaseq, seurat
---


# Question
Title: question to pick values based on conditions using a loop
I have a dataset with following information. Column C1\_1 is chromosome number, C1\_2 is SNP position and C3 is the p-value. I want to pick the most significant association within a genomic region of 3000 bp:

```
C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3

```

I want to split first column into two as:

```
C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L01 003000  L001g045    2.3e-6
L02 145000  L004g034    1.1e-2
L02 155000  L002g050    1.1e-4
L02 148000  L002g001    1.1e-3

```

Then I want to sort subsets such as `L01` and `L02` separately based on `C1_2` column:

```
C1_1    C1_2    C2  C3
L01 003000  L001g045    2.3e-6
L01 005000  L002g034    1.5e-12
L02 145000  L004g034    1.1e-2
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

```

Then I want to pick one value from each subset such as for `L02` – I want to pick the most significant value based on `C3` column under the difference of 3000. For example, `L02-148000` is chosen as shown below as its `C3` value is lower than `L02-145000`, and then the next value `L02-155000` is chosen as the difference is more than 3000 than the previous value of `L02-148000`:

```
C1_1    C1_2    C2  C3
L01 005000  L002g034    1.5e-12
L02 148000  L002g001    1.1e-3
L02 155000  L002g050    1.1e-4

```

I have a long dataset with 500,000 rows of values and the subset goes from `L01-L24`. I think it will need a loop to process these many values to get the final values.

This is just an example, I have shown one for L01 because the difference in L01 203000 and L01 205000 is less than 3000. I have shown two values for L02 because the difference between L02 145000 and L02 148000 is not more than 3000. Whereas the difference between L02 155000 and L02 148000 is more than 3000.

# Answer
> 4 votes
This question is deceptively challenging, but here is a potential solution:

```
library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_148000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_148000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 3 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   148000 L002g001 1.1e- 3
#> 3 L02   155000 L002g050 1.1e- 4

## additional example data
df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    2.3e-6
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 2.3e-06
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 5 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     5000 L002g034 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L02   159000 L002g050 1.1e- 4
#> 5 L03   165000 L002g001 1.1e- 3

```

Created on 2023-02-07 with reprex v2.0.2

## Edit

If you have two SNPs with the same C3 value (i.e. the top two rows in the example dataset below) you can select one of them using `slice_min(C3, with_ties = FALSE)`, e.g. the code below will select the 'first' entry when two C3 values are the same within the same group/subgroup:

```
library(tidyverse)

df <- read.table(text = "C1  C2  C3
L01_005000  L002g034    1.5e-12
L01_003000  L001g045    1.5e-12
L02_145000  L004g034    1.1e-2
L02_155000  L002g050    1.1e-4
L02_158000  L004g034    1.1e-2
L02_159000  L002g050    1.1e-4
L03_165000  L002g001    1.1e-3",
header = TRUE)

df
#>           C1       C2      C3
#> 1 L01_005000 L002g034 1.5e-12
#> 2 L01_003000 L001g045 1.5e-12
#> 3 L02_145000 L004g034 1.1e-02
#> 4 L02_155000 L002g050 1.1e-04
#> 5 L02_158000 L004g034 1.1e-02
#> 6 L02_159000 L002g050 1.1e-04
#> 7 L03_165000 L002g001 1.1e-03

df %>%
  separate(C1, into = c("C1_1", "C1_2")) %>%
  mutate(C1_2 = as.numeric(C1_2)) %>%
  arrange(C1_1, C1_2) %>%
  group_by(C1_1) %>%
  mutate(subgroup = cumsum(C1_2 > lag(C1_2, default = first(C1_2)) + 3000)) %>%
  group_by(C1_1, subgroup) %>%
  slice_min(C3, with_ties = FALSE) %>%
  ungroup() %>%
  select(-subgroup)
#> # A tibble: 4 × 4
#>   C1_1    C1_2 C2            C3
#>   <chr>  <dbl> <chr>      <dbl>
#> 1 L01     3000 L001g045 1.5e-12
#> 2 L02   145000 L004g034 1.1e- 2
#> 3 L02   155000 L002g050 1.1e- 4
#> 4 L03   165000 L002g001 1.1e- 3

```

Created on 2023-02-08 with reprex v2.0.2

# Answer
> 1 votes
### EDIT: Ignore this answer, it does not address one of your requirements - that of the 3000-long bin. You're going to need some partitioning operations which are not straightforward.

#### @jared\_marmot has cracked it in their solution.

This should be easy to do with dplyr/tidyr. Use

```
## Code not tested, might need tweaking
library(tidyr)
library(dplyr)
df_split <- df %>% separate(C1, sep = "_", into = c("C1_1","C1_2"))
df_aggr <- df_split %>% select(C1_1, C3) %>% group_by(C1_1) %>% summarise(C3 = min(C3)) %>% inner_join(df_split)

```

There are two problems that make your 3000 window complicated:

1. The start and end are not fixed: it's not 1 to n in steps of 3000; that'd be a straightforward bin.
2. We cannot even go with the smallest coordinate as the start, as you may have data points at 1500, 2500, 7200, 7600. If those data points exist, a set of 3000-long windows starting at 1500 would be 1500-4500, 4500-7500, etc. making your data go in 3 bins whereas your data points belong in 2 bins. You need windows with independent start coordinates. I've never encountered this requirement in my career even as a puzzle so I'm starting to think this might be an X-Y problem that needs to be addressed at a different level.

---
Tags: r
---


title: open and read Fasta file (raw data)

text:
# Question
Title: open and read Fasta file (raw data)
I have a big `fasta.dataset` file containing half a million proteins (1.0 GB). I have four lines for each protein code:

* line 1:the protein code
* line 2: protein length in amino acids
* line 3: amino acid sequence
* line 4: secondary structure

Now, I am trying to open and read it in python (Biopython), and it does not work:

```
filename = 'pdb.fasta_qual.dataset'
sequences = SeqIO.parse ( filename,'fasta')
for record in sequences:
    example = record
    break
print(example)

```

How can I read it in python and loop through the file to look at line 3 for each protein to count the sequence length and distribution?

here is the first 5 lines of my file: which my file contains 500,000 proteins for each one has a 4 lines (name ,len of protein in amino acid,the seq represents by letters which what I would to calculate,the secondary structure)

```
4LGTD
247
M       S       E       K       L       Q       K       V       L       A       R       A       G       H       G      T
.       .       E       E       H       H       H       H       H       H       H       T       T       S       S       .

```

I want to open and read the file and loop through line 3 for each protein to calculate the length of the sequences and plot a histogram ! to check the distribution. ,,,

The output i am expecting is :

The len for the first seq is =

The len for the second seq is =

Until

the len of the last sequence which is number (500.000)=

And then i can plot a histogram for the len of the sequences ,,,

NOTE: I have opened and read the file's info by Linux, but I could not by python.

# Answer
> 0 votes
I better try and explain the issues a bit better. A dataset of 0.5 million at 1 GB will struggle to load into `Biopython`. `Biopython` isn't great at handling large numbers of sequences - for a regular desktop/laptop - and what will happen is the parser will freeze or else take an enormous amount of time (many hours). I call this a "RAM bottleneck". A generator is an easy workaround.

The other issue with the code is there are lots of bugs combining format and code. There is a clarity issue as well it looks like you want to calculate the secondary structure - rather than import the secondary structure alongside the alignment. I could be wrong on this though, certainly secondary structure is present in the *ad hoc* alignment. There is an alignment format that handles both secondary structure and sequence data. You might think about restructuring the question and your approach into multiple questions.

1. Reformatting the *ad hoc* alignment to permit import into `Biopython` on a *pilot data set* (not 0.5 million sequences). This is close to @terdon's answer.
2. Code that will process very large files in `Biopython`. This is the answer below.
3. An alignment format via `Biopython` (probably output - could be input) that integrates linear sequence with secondary structure.
4. Plotting a histogram for 0.5 million sequences.

Points (Questions) 1. to 3. can be leveraged by Biopython. Point (Question) 4 can be done via Python too - graphing/histograms are pretty good. There seems an underlying issue of project management for coding based work.

If this is a RAM bottleneck then solution is below.

---

The easiest approach is to use a generator (please see 'format').

```
example = (record for record in SeqIO.parse(filename,'fasta'))
print (next(example)) # example.id ?

```

A generator will deal with the sequences one by one instead of trying to load the whole thing to memory. It may not be suitable for certain complex calculations, but it is by far the easiest approach to the problem (if the problem is the OP is struggling to load a million sequences into memory).

**Note** I always used `yield` commands but @Steve introduced the `()` format and I definitely prefer it. Also note that manipulating the output of a generator can be different to a list, tuple or dictionary (calculation depending).

---

**Format** The format supplied isn't fasta so the `SeqIO` will not work in `fasta` mode and should throw an error/exception. The format you need to leverage `Biopython` is very different. I would recommend a standard 'open' command and passing this to a custom dictionary. Also the output of `print(next(record))` should produce an object symbol, not the actual data e.g. `record.id` would be needed to see data: it does however demonstrate the generator is working.

# Answer
> 4 votes
It sounds like all you really want is to print every 2nd line of a text file. If so, you don't need Python, let alone BioPython, you can do it with basic \*nix tools:

```
$ awk 'NR%4==2' pdb.fasta_qual.dataset
247

```

`NR` is the current line number, and `%` is the modulo operator. Therefore, `NR % 4` will equal `2` on the second line of each group of 4 lines in the file. In `awk`, the default action when something evaluates to true is to print the current line, so this will print out all length lines from your input file (assuming you don't have blank lines anywhere, which we don't know since you only showed the first entry).

If you want to add some sort of message to be printed, you could do:

```
$ awk 'NR%4==2{ print "The length of protein",++c,"is:",$0}' pdb.fasta_qual.dataset
The length of protein 1 is: 247

```

If you really must do this in Python, you could do something like this:

```
#!/usr/bin/env python3

import sys

#my_file= sys.argv[1]
line_number = 0
with open(sys.argv[1]) as f:
    for line in f:
        line_number += 1
        if (line_number % 4 == 2):
            print("The length of protein %d is: %d" % (line_number,int(line.strip())))

```

And then you would run it like this (assuming you named the script `foo.py`):

```
$ python3 foo.py pdb.fasta_qual.dataset
The length of protein 2 is: 247

```

---

Out of curiosity, I created a test file by repeating your example 10000000 times, which left me with a 2.4G file:

```
$ perl -e 'open(A,"pdb.fasta_qual.dataset"); @lines=<A>; print "@lines" x 10000000' > testFile
$ ls -lh testFile
-rw-r--r-- 1 terdon terdon 2.4G Feb  7 17:06 testFile

```

I then tested both approaches for speed. Unsurprisingly, `awk` was much faster:

```
$ time awk 'NR%4==2' testFile > /dev/null

real    0m6.447s
user    0m5.917s
sys     0m0.492s

$ time foo.py testFile > /dev/null

real    0m10.856s
user    0m10.111s
sys     0m0.523s

```

---
Tags: fasta, proteins, protein-structure, biopython
---


# Question
Title: How do I change the legend for a Violin Plot with median dot
I'm plotting scRNA data using Violin Plots. I need to re-order the plots and add median dots. My code is:

```
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9))

Vln_plot$data$split<- factor(x = Vln_plot$data$split,
                                 levels = c("WT", "F7KO", "F8HET"))

Vln_plot

median.stat <- function(x){
  out <- quantile(x, probs = c(0.5))
  names(out) <- c("ymed")
  return(out)
}

```

and my output looks like:

Is there any way to fix the legend so those dots aren't in each box? I also need to figure out how to get rid of that x-axis tic mark and change the x axis label from Identity to something else. Can anyone point me in the right direction on how to accomplish this?

# Answer
> 1 votes
Your command looks like it creates a ggplot2 object, so it should be possible to remove the symbol via `guides`, setting the associated aesthetic to `"none"`:

```
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9)) +
  guides(color = "none")

```

If that doesn't work, try telling the `stat_summary` command to disable the legend via `show.legend = FALSE`:

```
Vln_plot <- VlnPlot(object,
                    features = c("Mcidas"),
                    cols = c("coral1", "cadetblue", "cadetblue2"),
                    split.by = "genotype",
                    assay = "RNA",
                    pt.size = 0,
                    log = TRUE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9), show.legend = FALSE)

```

# Answer
> 1 votes
If you can provide a minimal reproducible example with a 'cut-down' version of your data it would be a lot easier to answer your question. I have an example from a similar question on stackoverflow (https://stackoverflow.com/a/73696447/12957340) that illustrates one approach to your problem:

```
# Load libraries
library(tidyverse)

# Load example data
raw_data <- structure(list(Type = c("EY", "EY", "EY", "EY", "EY", "EY", "EY",
                                    "EY", "EY", "EY", "EY", "EY", "EO", "EO", "EO", "EO", "EO", "EO",
                                    "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EY", "EY",
                                    "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EY", "EO",
                                    "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO", "EO",
                                    "EO", "EO", "EO", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe",
                                    "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe", "Fe"),
                           Maturity = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L,
                                                  1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
                                                  2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Before",
                                                                                "After"), class = "factor"), Mn = c(49.65414679, 30.42267647,
                                                                                53.57468125, 58.38531743, 11.07088615, 45.10644112, 41.33779076,
                                                                                18.77426257, 60.36149117, 67.0260032, 38.13116392, 69.58713059,
                                                                                66.80825471, 30.34496525, 56.90875291, 72.37427455, 22.06855022,
                                                                                39.26030375, 35.31067394, 24.95698024, 24.88586385, 34.65821779,
                                                                                56.2750941, 19.56655377, 45.04117517, 44.04604616, 30.4130276,
                                                                                15.76824315, 21.50636001, 27.25777212, 17.08214603, 14.17493713,
                                                                                34.24593956, 12.39548126, 23.59787998, 38.42021152, 32.54843776,
                                                                                18.47119505, 26.05687922, 32.97977895, 16.91703485, 30.0626039,
                                                                                26.72724955, 17.69433512, 22.34314521, 16.73815161, 13.08749188,
                                                                                19.60771131, 24.08865106, 34.57472809, 14.94650052, 32.60745327,
                                                                                29.86439635, 20.71564742, 22.32568148, 18.37220522, 24.22830773,
                                                                                20.91660251, 20.96486868, 34.25266591, 35.57200211, 14.5736099,
                                                                                25.67187037, 27.7352539, 42.919751, 27.1949762, 18.94238981,
                                                                                25.25867502, 31.66215867, 18.1734287, 24.41910112, 12.07438938,
                                                                                12.85034127, 19.75012923, 35.96030792, 30.07986478, 38.91460672,
                                                                                41.18492241, 41.68085289, 60.93940724, 49.77930187, 17.56806894,
                                                                                37.13306953, 47.06926881, 85.62093648, 24.31500546, 22.049633,
                                                                                40.42923285, 87.53339143, 37.08699364, 38.29525586, 28.00272281,
                                                                                15.58855204, 37.14164977)), row.names = c(NA, -94L), class = c("tbl_df",
                                                                                "tbl", "data.frame"))

```

---

Functions:

```
summary_df <- raw_data %>%
  mutate(Maturity = ifelse(Maturity == "Before", "a", "x")) %>%
  group_by(Type, Maturity) %>%
  summarise(max_val = max(density(Mn)[[1]]) + 4) %>%
  ungroup()

median.stat <- function(x){
  out <- quantile(x, probs = c(0.5))
  names(out) <- c("ymed")
  return(out)
}

```

---

Similar 'problematic' violin plot:

```
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA, position=position_dodge(.9)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9))
#> Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
#> ℹ Please use the `fun` argument instead.

```

---

Remove dots:

```
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA, position=position_dodge(.9)) +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(fun.y = median.stat, geom='point', size = 10, colour = "blue",
               position = position_dodge(0.9), key_glyph = "blank")

```

---

Remove dots, but keep 'colour' legend to specify what the dots represent:

```
ggplot(raw_data, aes(x = Type, y = Mn, fill = Maturity)) +
  geom_violin(trim = FALSE, position="dodge") +
  geom_boxplot(width=.1, outlier.shape=NA,  position=position_dodge(.9), key_glyph = "boxplot") +
  scale_y_continuous(breaks = seq(0, 100, 25)) +
  theme_classic() +
  theme(axis.title.x = element_blank(),
        axis.ticks.x = element_blank()) +
  labs(y="Percent Maturity", x= "") +
  geom_text(data = summary_df,
            aes(label = Maturity,
                y = max_val,
                x = Type,
                group = Maturity),
            position = position_dodge(width = 0.9),
            check_overlap = TRUE,
            show.legend = FALSE,
            inherit.aes = FALSE) +
  stat_summary(aes(colour = "Median"), fun.y = median.stat,
               geom='point', size = 10,
               position = position_dodge(0.9)) +
  scale_color_manual(values = "blue", name = "Statistics") +
  guides(fill = guide_legend(override.aes = list(shape = NA)))

```

Created on 2023-02-08 with reprex v2.0.2

Hopefully this approach can be adapted to your use-case.

---
Tags: r, seurat
---


# Question
Title: how to do manual search on NCBI with biopython and add link to publication
I'm trying to reproduce a complex manual search on NCBI using BioPython, but obviously something is not right. My command line is not working, but I'd especially like a way to have an accurate search where I can select only the genomes related to a published paper. Is there an easy way to do that?

I want to limit my results to only those genomes that i) have an NCBI BioProject entry and ii) that BioProject entry contains a `Publications` field.

For example, let's say I want all the genomes of Microbacterium, the latest version, from RefSeq and not partial (this is what I tried to write but something is not working) and I would like to add one more filter. Take, among the remaining/selected genomes, only the ones related to published papers.

For this strain, we have a related paper: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA669461/

For this one, we do not: https://www.ncbi.nlm.nih.gov/bioproject/PRJNA608230/

and I really would like to not have to check 100 genomes manually, but I don't know how to manage it.

I tried:

```
query = "Microbacterium[Organism] AND latest refseq[filter] AND (all[filter] NOT
partial[filter])"
handle = Entrez.esearch(term=query, db="Assembly", retmax=600)
ids = Entrez.read(handle)["IdList"]

```

And after this, I'm just using a `for` loop to get genomes, but my query is not good and it is not handling the problem of the publication.

Thanks for any help!

---

Edit : as I need to search in the db Assembly unfortunately the answer with the flag author or date is not working

# Answer
Entrez has the tags \[journal\], \[volume\], \[page number\], \[title\], \[year\] and \[author\] to identify all sequences from a given publication. You will find a full list of tags here from the drop down menu.

Thus,

```
query = "Microbacterium[Organism] AND Nature[Journal] AND 528[volume] AND Bai[author] AND 364[page number]"

```

This will identify all *Microbacterium* genome sequences associated with Bai *et al* (2015) Nature 528:364- and is here.

For practical purposes I routinely search using the organism and author tags for their sequence output.

---

The example given (post-edit) was

> Kim KK et al., "Microbacterium aquimaris sp. nov., isolated from seawater.", Int J Syst Evol Microbiol, 2008 Jul;58(Pt 7):1616-20

```
query = "Microbacterium[Organism] AND 58[volume] AND Kim[author] AND 1616[page number] AND 2008[year]"

```

I removed the \[journal\] tag because the abbreviation it uses is a bit unusual and replaced it with \[year\]. Still works here.

---

From the comments ... if only published sequences were required then using the `range` function and `[publication date]` is one way in ...

```
query = "Microbacterium[Organism] AND 2000/1/1:2020/1/1[publication date]"

```

The best strategy is to use the `[author]` because data production tends to be dominated by a small number of groups.

---

**Comments of Comments**

Why isn't searching working. The following search works

```
query = "Microbacterium[Organism] AND (all[filter] NOT partial[filter])"

```

The term it doesn't like is `latest refseq[filter]` if that is removed it will work. `refseq` works, but it doesn't like `latest`.
> 0 votes

---
Tags: python, biopython, ncbi
---


# Question
Title: Arlequin - How to assign samples to populations?
I have converted a VCF file to the arlequin file format. It consists of 20 samples that *should* belong to three populations, along with a few thousands SNPs. However, when I load the arlequin file into the software, it shows I only have 1 population. Any way to fix this?

Thanks.

# Answer
Arlequin format is historically extremely difficult to manually manage. The group leader eventually wrote software to allow easy Arlequin format management and conversion. This is achieve by Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018). All information is available here and the manual is here. I would personally avoid attempting to manually rectify population designation.
> 0 votes

---
Tags: snp, format-conversion
---


title: Nextflow tuple from CSV cols

text:
# Question
Title: Nextflow tuple from CSV cols
I am trying to create a channel with a CSV file and use it as tuple input. Yet the output I obtain does not contain the entire file but only the first element of the file (1st row).

**the csv**

```
dig,val
1,One
2,Two
3,Three
4,Four
5,Five

```

output : `cat digOne.txt` `1`  instead of the whole collection digOne, digTwo...

```

 size_amp_ch = Channel
        .fromPath( params.file)
        .splitCsv(header:true).map{
            row  -> tuple(row.dig,row.val)
        }

process{
    input:

    tuple val(digits), val(values)

    output:

    path "dig_${values}.txt"

    script:

    """
    echo $digits >> dig_${values}.txt

}

```

thank you

# Answer
> 1 votes
I think the problem is just that you have extra whitespace at the end of the first line in your CSV file. Either remove the whitespace manually, or use the `strip` parameter with the `splitCsv` operator to remove leading and trailing whitespace from values. For example:

```
params.file = './file.csv'

process test_proc {

    input:
    tuple val(digits), val(values)

    output:
    path "dig_${values}.txt"

    """
    echo "${digits}" > "dig_${values}.txt"
    """
}

workflow {

    Channel.fromPath( params.file ) \
        | splitCsv( header:true, strip:true ) \
        | map { row ->
            tuple( row.dig, row.val )
        } \
        | test_proc \
        | view()

}

```

Results:

```
$ nextflow run main.nf
N E X T F L O W  ~  version 22.10.0
Launching `main.nf` [stupefied_boyd] DSL2 - revision: 2a9d8dc414
executor >  local (5)
[4f/a73196] process > test_proc (2) [100%] 5 of 5 ✔
/path/to/work/9d/4905d0236788c1a3733da5358b7fc4/dig_Four.txt
/path/to/work/f9/d04e8a1c2265d73cde626703317df7/dig_Five.txt
/path/to/work/ff/522238f87a5641fdc7adce3e4d1cea/dig_Three.txt
/path/to/work/e9/9b7ec2a65d45986fd8944da4a2054c/dig_One.txt
/path/to/work/4f/a73196fbd56db3c9647fa7b072edfd/dig_Two.txt

```

---
Tags: nextflow, workflow-management
---


# Question
Title: How much does Nanopore cDNA Sequencing Cost?
*\[this question is based on a question that was asked on Reddit\]*

We're interested in doing whole-transcriptome cDNA sequencing of 24 mouse cell lines, and are deciding between Illumina and Nanopore sequencing. What are the costs involved in Nanopore sequencing, and how much hassle is it to analyse Nanopore data?

# Answer
> 6 votes
The standard PCR-cDNA barcoding kit is \\$750 USD, and has enough reagents for six runs, with each run using up to 24 samples (i.e. 144 samples total). Each run also needs one MinION flow cell (\\$900 USD).

The total cost including external reagents is naturally dependent on the number of samples, and gets cheaper as the number of samples increases. For 24 samples, it's about \\$5200 NZD (\\$3100 USD). See below for more details.

If you could outsource this to a lab that has a PromethION (or P2 Solo), the sequencing would be a great lot cheaper. PromethION flow cells start at \\$3600 for a pack of 4 (i.e. the same unit cost as MinION flow cells), but have about ten times the yield of MinION flow cells (which means more samples, or greater depth).

From a bioinformatics perspective, there's more work to do with nanopore reads because there are no commonly used / established pipelines (so I've had to create my own), but the results are much cleaner once they've been processed.

**Cost Breakdown**

For cDNA sequencing in Aotearoa, MinION works out cheaper than Illumina when compared against a service centre, about two thirds the price for us. Here is the reagent list I wrote up a couple of months ago for sequencing from 30 whole transcriptome samples, mouse or human (prices in NZD). This is with an ideal read threshold of 1M reads per sample (which I have found to give good quantitative results for differential expression). You could probably get away with lower yield for microbial sequencing:

ONT reagents:

* 1 x PCR-cDNA Barcoding Kit 24 \\$750 USD \[\\$1210 NZD\]
* 1 x Rapid Adapter Auxiliary V14 \\$288 USD \[\\$481 NZD\]
* 1 x R10.4.1 flow cell @ \\$900 USD \[\\$1503 NZD\]

Total ONT cost: \\$3194

Other reagents:

* 10 mM dNTP NEB N0447 \[\\$145 NZD\]
* Maxima H Minus RT ThermoFisher EP0751 \[\\$346 NZD\]
* RNAclean or AmPure XP beads \[AmPure XP: \\$800 NZD\]
* Exonuclease 1 NEB M0293 \[\\$160 NZD\]
* RNAseOUT ThermoFisher 10777019 \[\\$569 NZD\]

Total external reagent cost (excludes nuclease-free water, tubes, pipette tips, gloves, ethanol, DNA QC): \\$2020 NZD

\[As well as whatever you need to extract 200ng total RNA per sample\]

Our current protocols for this PCR-cDNA barcoding + V14 RA approach can be found here.

The comparative Illumina cost according to the service centre we asked a couple of months ago was \\$185 NZD per sample + \\$5000 NZD per flow cell (or $9440 NZDfor 24 samples).

For almost all of these external reagents the amount ordered is an absolute minimum (from memory I think most are something like 50 reactions), so doing a second equivalent experiment would mostly only need the flow cells to be purchased.

In addition to the cost, turnaround time with MinION is quicker (for us it's a couple of weeks versus a couple of months for Illumina \[with most of that time spent waiting in a queue\]), and the results are nicer to analyse: full length transcript sequences with very even coverage profiles.

I'm aware that sequencing is more expensive in Aotearoa than it is in other countries, so the cost / time difference may not be as extreme elsewhere.

**Bulk discounts for ONT flow cells**

The prices above are for ordering reagents and flow cells one at a time, using the most expensive reagents. I call this a "grant funding" budget, because it's the quote I would give to people who are looking to fund research. There are ways to reduce costs when organising sequencing in bulk, and I've already hinted above at the extra capacity in reagent kits.

However, the biggest savings for nanopore sequencing can be achieved by ordering lots of flow cells at once. For example, ordering 48 MinION flow cells means that the per-flowcell price is \\$500 USD instead of \\$900 USD. That's 48 flow cells spread over the course of an entire year, ***with costs paid on delivery***. If you can encourage multiple research groups to put in a combined bulk order, it's not too far out of the realm of possibility.

**The long-read advantage**

I've been told that a researcher is able to get 50M paired-end Illumina reads from a sequencing service centre for about $100 CAD.

50M paired end reads is about 10 Gb, which can be achieved on a MinION run (\\$600 - \\$1000), or about 1/10th of a PromethION run (\\$850 - \\$1500).

Even given those prices, it doesn't take much of a long-read advantage for MinION to become the most economical option. My own studies suggest that 1M full-length reads from nanopore sequencing (i.e. about 1/10th of a MinION R9.4.1 flow cell) gives results that are at least as good as 20M short paired-end Illumina reads.

---
Tags: rna-seq, nanopore, illumina
---


# Question
Title: How do I measure the proportion of different bacteria in a sample from a high-throughput sequencer?
*\[this question is based on a question that was asked on Reddit\]*

I have sequenced some \[mostly cell-free\] DNA from a sputum sample on a nanopore sequencing machine, and would like to know what is in it (i.e. the proportion of different bacteria). The run produced 2 million reads with an average read length of about 200bp.

What's the best way to go about doing a metagenomic classification of sequenced reads?

# Answer
> 4 votes
If you've got a fast solid state drive with \>200G of space, then kraken2 + bracken works really well on both long and short reads (e.g. Nanopore, Illumina). It's also specifically designed for metagenomic analysis, unlike other general aligners like STAR, bwa, bowtie2, etc., so you get tree-based classification which can deal with partial / incomplete matches.

First, use kraken2 to map reads to taxa in memory-mapped mode (to reduce system memory consumption):

$ kraken2 --threads 10 --db /mnt/ultra\_fast/kraken/metagenome --quick --memory-mapping reads/\*.fastq.gz --report kraken2.report.txt \> kraken2.out.txt

Then, use bracken to adjust the counts based on mapping / library probabilities:

$ ~/install/bracken/bracken -d /mnt/ultra\_fast/kraken/metagenome -i kraken2.report.txt -o bracken.report.txt

I use pavian to look at the results, because it creates quite understandable Sankey plots:

But inspecting the `_bracken_species.txt` file is a pretty good text-based way to do it. If there's substantial amounts of a particular bacteria in the sample, it should come out in the results:

```
100.00  232939  0       R       1       root
99.99   232912  0       R1      131567    cellular organisms
68.52   159618  0       D       2           Bacteria
67.67   157635  0       P       1224          Proteobacteria
67.47   157152  0       C       1236            Gammaproteobacteria
61.82   144001  0       O       91347             Enterobacterales
61.79   143941  0       F       543                 Enterobacteriaceae
32.28   75183   0       G       561                   Escherichia
32.26   75147   75147   S       562                     Escherichia coli
0.02    35      35      S       1499973                 Escherichia marmotae
0.00    1       1       S       564                     Escherichia fergusonii

```

kraken2/bracken indexes can be downloaded from here (the build process takes *ages*; it makes more sense to use the pre-built indexes):

https://benlangmead.github.io/aws-indexes/k2

As far as I'm aware, kraken2 can't do strain-level assignment with the most commonly-used databases (although that may depend on what information that is stored in the database). For strain-level assignment, it seems like MetaMaps may work:

https://doi.org/10.1038/s41467-019-10934-2

---
Tags: sequence-analysis, nanopore, metagenome, illumina
---


# Question
Title: Is there a way to calculate the number of parsimony informative characters in a morphological dataset in R?
I have a dataset of binary morphological characters and would like to know how many of them are parsimony informative. Is there a way to calculate this with an R package?

# Answer
The `ape` library (Analyses of Phylogenetics and Evolution) using the `MPR` method (Most Parsimonious Reconstruction).

Information on `ape` is here.

The format

```
install.packages(“ape”)
library(ape)
MPR(x, phy, outgroup)

```

> **x** : a vector of integers.

> **phy** : an object of class "phylo"; the tree must be unrooted and fully dichotomous.

> **outgroup** : an integer or a character string giving the tip of phy used as outgroup.

Rooting the tree is recommended. In the example rooting isn't used but an outgroup is used as taxa `t1`, that will be fine.

Generating random data,

```
vectx <- rpois(50, 1) # make random data
tree <- rtree(50, rooted = FALSE)  # also rtree(50) make random tree
MPR(vectx, tree, "t1")

```

Thats it, there's nothing more. Data, tree = steps.

**Clarification** the question is number of parsimonious characters. As a strict definition the number of informative characters is the same as the number of steps and thats what `MPR` will describe ... thats easy.

There are other possible meanings

1. The number of parsimonious characters/steps to make a tree is calculated via `MPR` as above.
2. if its the actual changes at a given position in the alignment (or vector) for each node in the tree - thats the answer below - but thats complicated.
3. If its simply is this site (not character) is parsimoniously informative and this one isn't then the "alignment" or vector could be used one character position at a time and subject to `MPR` and the step length would determine how parsimoniously informative each site was ... thus point 3 is easily doable in `R`

---

**Note** Calculating the actual sites on the alignment and associated changes per node is easily done outside `R`, but within `R` I don't know if thats doable.

Wait, I'm wrong it's doable but I would envisage its difficult to implement. Calculating parsimonious characters for each given node in a tree - yes I do that - doing that in `R` ... errhhh thats not me. There are Pythonic ways in, but I've not implemented them.
> 1 votes

---
Tags: r, phylogenetics
---


# Question
Title: calculate the backchain RMSD between two pdb files by pymol in Python
i have two protein pdb files and want to calculate the backchain RMSD between them. As far as I know the GUI of Pymol can use align to calculate the RMSD (but I don't know if it is the backchain RMSD).I want to use the python code to use pymol to calculate the backchain RMSD. I know how to align them but don't know how to calculate the RMSD.

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

```

# Answer
The backbone can be specified using a selector for the names of these, which are `N`, `CA` (Cα) and `C` (C-prime).

```
out = pymol.cmd.align('foo and name N+CA+C and polymer', 'bar and name N+CA+C and polymer')
rmsd, n_atoms, n_cyles, n_rmsd_pre, n_atom_pre, score, n_res = out

```

If you do not wish the proteins to be moved you can use `rms_cur`, which is just the rmsd.
> 2 votes

---
Tags: python, pymol
---


title: Energy minimization of molecules

text:
# Question
Title: Energy minimization of molecules
I have 3000 molecules which I need for docking but before that I need to perform energy minimization for all the molecules and it is taking a long time because I am doing it one by one. So is there a way that I can select multiple molecules and do the energy minimization

# Answer
> 1 votes
As @M\__ has mentioned, the terminology in your Q is a bit confusing.

A molecule can be the small molecule (ligand, compound etc.) or a polypeptide chain (protein).

Conformer generation of a compound is commonly not called 'energy minimisation' for technical reasons. For example, you normally generate in rdkit a conformer via `AllChem.EmbedMolecule(mol)`, but you could do it via

```
p = AllChem.MMFFGetMoleculeProperties(mol, 'MMFF94')
ff = AllChem.MMFFGetMoleculeForceField(mol, p)
m = ff.Minimize()
assert m == 0, 'failure'

```

The latter would make sense if you want to add constraints (e.g. `ff.MMFFAddPositionConstraint(0, 0.1, 10)`).

There is also DFT (e.g. Gaussian or Schrödinger Jaguar), wherein the conformer is generated from a forcefield, but using a slimmed-down quantum theory to described orbitals as opposed to atoms. This is much slower than the above which use atoms as rigid particles tethered together (molecular mechanics).

The `EmbedMolecule` uses ETKDG method, which is a stochastic search method that cheats by using knowledge derived from crystal structures (cf. CCDC), so is way faster. It will be based on crystal structures, where the latter makes the compound best suit the terms of the MMFF94 forcefield.

If you are using a conformer generator, or similar, then the only quicker method is to retrieve the crystallographical structure.

If you are talking about minimising a whole protein, every forcefield will give slightly different results and there is no DB you could use. However, if it is for a snapshot/pose calculation, a common hack is do a local minimisation around your target and use implicit solvent which will be much faster. You have added the tag `MD`, which would suggest FEP calculations, for this local minimisation would not work.

---
Tags: docking, molecular-dynamics
---


# Question
Title: Get a comprehensive list of virulence factors for certain or all pathogens?
Is there a database with all known proteins/genes of specific pathogens such as Staph aureus and annotated whether it is a known virulence or resistance factor that idealy can be downloaded in table format?

# Answer
> 4 votes
I suggest that you are looking for a resource like VFDB (virulence factor database). That table does annotate by bacterial species, so you can filter as you like.

You can download a description table or FASTA sequences from their downloads page (see screenshot). Last update was 4 days ago, so they should be quite current.

For resistance genes, you could use something like CARD (Comprehensive Antibiotic Resistance Database), which also has a downloads page where you can download a similar TSV table or FASTA sequences.

Both resources include web tools for predicting the presence of the elements of interest in sequences.

If you are looking for a state of the art prediction tool for such elements in genomes, you might look into seqscreen from Todd Treangen's lab.

For a general purpose microbiological annotation and analysis resource, you might also look into BV-BRG, which succeeds PATRIC. I have only used PATRIC but BV-BRG seems very thorough.

I think that you are unlikely to find a resource that has *non-virulent/non-resistance* genes annotated.

---
Tags: public-databases, microbial-genomics, pathogenicity
---


title: Primer trimming-fasta files

text:
# Question
Title: Primer trimming-fasta files
My goal: I want to trim off the primers (Forward : CGAGAAGACCCTRTGRAGCT, Reverse : GTTGGGGYGACCNYGG) from a fasta file with a lot of dna sequences allowing for some (e.g. 3) mismatches (identity).

Input sequences Seq1:

> **CGAGAAGACCCTATGGAGCT**TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG**GTTGGGGCGACCGCGG**AACGT

Seq2:

> ATGGCCCATCCCTCACAGCTAGGATTTCAAGATGCAGCTTCCCCAGTTATAGAAGAACTTCTCCACTTTCACGACCATGCCCTAATAATCGTTTTTTTAATTAGTACACTAGTA

Seq3:

> CATAAGA**CGAAAAGACCCTATGGAGCT**TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG**GTTCGGGCGACCGCGG**

Output sequences

Seq1:

> TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG

Seq3:

> TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

Seq2 is absent because it has no match for the primers so it should be discarded.

In seq 3 there is a mismatch in Forward (CGA**A**AA...) and in Reverse (GTT**C**GG...). I don't need the primers to have 100% identity. 85-90% is OK.

Issues:

1. primer sequences can be anywhere in each sequence, not just in the beginning.
2. The sequence between the primers is not constant. it is approximately 80 bases but not exactly 80.
3. I don't know how to put the primers' identity parameter.
4. How can I search for R or N nucleotides (IUPAC)? Because I don't know how, I wrote the primers with specific nucleotides and not R or N.

This is the code I have so far, but it has the problems mentioned above:

```
from Bio import SeqIO

def trim_adaptors(records, adaptorF,adaptorR):

    len_adaptorF = len(adaptorF)
    len_adaptorR = len(adaptorR)

    for record in records:
        indexF = record.seq.find(adaptorF)
        indexR = record.seq.find(adaptorR)
        if indexR == -1  or indexF == -1:
            # adaptor not found, so won't trim
            print(record.id, "no primer/s match")
            continue
        else:
            # trim off the adaptor
            yield record[indexF + len_adaptorF:indexR + 1]

original_reads = SeqIO.parse("trimming_testfas.fas", "fasta")

trimmed_reads = trim_adaptors(original_reads,"CGAGAAGACCCTATGGAGCT" ,"GTTGGGGCGACCGCGG")

count = SeqIO.write(trimmed_reads, "trimmed.fasta", "fasta")

print("Saved %i reads" % count)

```

This doesn't have to be a script in Python, I am also open to using existing tools.

Thank you a lot in advance, hope you got me :)

# Answer
The IUPAC ambiguity codes can be thought of as regular expression character classes. `R` matches any purine, so `[AG]`, `Y` matches any pyrimidine, so `[CT]`, and `N` matches anything at all, so `.` or `[ACTG]`. So those are easy to handle: instead of  `CGAGAAGACCCTRTGRAGCT` you can use the regular expression `CGAGAAGACCCT[AG]TG[AG]AGCT`, and instead of `GTTGGGGYGACCNYGG` you can use `GTTGGGG[CT]GACC.[CT]GG`.

Allowing for random mismatches is trickier though and since you're open to using existing tools, I strongly urge you to use something like `cutadapt` instead of trying to reinvent the wheel. So, install `cutadapt` on your system (see here for instructions) and then run:

```
cutadapt -a "CGAGAAGACCCTRTGRAGCT...GTTGGGGYGACCNYGG;max_error_rate=0.15;" --discard-untrimmed -o output.fa file.fa

```

That assumes your input sequences are in fasta format in the file `file.fa` and will create the output file `output.fa`:

```
$ cat output.fa
>seq1
TAAGGCGCCAGAACAGCTCACGTCAAACACCCCCGCATAAAGGGAATAAACCAAGTGGACCCTGCTCTAGTGTCTTTG
>seq3
TTAGACGTCAGAGCAGCTCATGTAAAGCACCCCTAAACAAAGGAAAAAACCAAATGAAATCTGCCCTAATGTCTTTG

```

The options used are `-a $primer1..$primer2` which tells `cutudapt` these primers are linked adapters (https://cutadapt.readthedocs.io/en/stable/guide.html#linked-adapters) and then we also add `max_error_rate=0.15;` to allow up to 85% identity matches. Cutadapt already knows how to handle IUPAC ambiguity codes so nothing special is needed for those.

I am not sure how this will handle cases where only one of the two sequences is present, you'll need to test for that. I expect that they will be discarded since we're specifying linked primers but I would double check to be sure. In any case, `cutadapt` is a very powerful tool, so I'm sure there will be a way to handle that in the way you want if you check the documentation.
> 2 votes

---
Tags: python, trimming, primer, multi-fasta
---