# Question
Title: Maintaining SNP names when calculating allele frequencies from a VCF file on plink 1.9
I downloaded the following vcf file from https://www.internationalgenome.org/.

http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chr22.phase3\_shapeit2\_mvncall\_integrated\_v5b.20130502.genotypes.vcf.gz

I wanted to first try and calculate allele frequencies, and then from there calculate allele frequencies based on populations and so on. However, I ran into an issue where while I was using plink, I would utilize the --vcf command to try and calculate frequencies. The frequencies were calculated, but the SNP names were missing. I used the following command in Git Bash on Windows:

```
    plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.gz --freq --out chr22datanew

```

I then tried to convert the file to plink binary (bed, bim, fam) to try and use that for the basis of the allele frequencies. I used the following commands:

```
plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes.vcf.vcf --make-bed --out chr22data

plink --bfile chr22data --freq --out chr22freqs

```

Once again, the frequencies were calculated with no issue, but the SNP names were missing. I checked the .fam file that was created, and the SNPs were there, but not in the .bim file. This may be a simple, elementary level question, but does anyone know how to fix this? I am using plink 1.9 instead of vcftools because I'm using a computer that runs on Windows OS.

# Answer
> 1 votes
The standard and most common way to name SNPs is by using rsIDs.

The reason your output from PLINK does not include SNP names is because these are not included in the original VCF you are working with. There are at least two possible options in order to fix this.

1. Annotate first the source VCF with rsIDs from dbSNP. You can do this using bcftools and a recent (or latest) human rsID compilation in VCF format. E.g.

> Downloading rsIDs

```
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/All_20180423.vcf.gz
wget ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/All_20180423.vcf.gz.tbi

```

> Subsetting chr22 from annotations

This makes the annotation of the chr22 VCF faster

```
bcftools view --threads 10 --regions 22 All_20180423.vcf.gz -Oz -o 22_20180423.vcf.gz
bcftools index 22_20180423.vcf.gz

```

> VCF Annotation

```
bcftools annotate --threads 10 -a 22_20180423.vcf.gz -c ID -o ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes_annotated.vcf.gz

```

> PLINK frequency estimations

Now you can proceed as you had done already

```
plink --vcf ALL.chr22.phase3_shapeit2_mvncall_integrated_v5b.20130502.genotypes_annotated.vcf.gz --make-bed --out chr22data
plink --bfile chr22data --freq --out chr22freqs

```
2. An alternative option is to build your own 'SNP names' based on chromosome and cordinate (e.g. 22\_16255382) using plink commands.

---
Tags: vcf, snp, plink, allele-frequency, rsid
---


# Question
Title: Alignment process in Nextflow not picking up the index files from the index process
I'm using bowtie2 to build a small index. After building the index, I want to pass it to the bowtie2 alignment process. The issue is bowtie2-build outputs multiple index files with the .bt2 extension. When I try to provide these index files as an input to the Alignmnet process, it gives the following error:

> Error executing process \> 'SAM\_FOR\_STRAND (1)'
>
> Caused by: Process `SAM_FOR_STRAND (1)` terminated with an error exit status (255)
>
> Command executed:
>
> bowtie2 -t -x index.1.bt2 index.2.bt2 index.3.bt2 index.4.bt2 index.rev.1.bt2 index.rev.2.bt2 -1 SRR5655560\_1.atria.fastq -2 SRR5655560\_2.atria.fastq -S "SRR5655560.sam"
>
> Command exit status: 255
>
> Command output: (empty)
>
> Command error: (ERR): "index.1.bt2" does not exist or is not a Bowtie 2 index Exiting now ...
>
> executor \> local (4) \[71/c095e4\] process \> INDEX \[100%\] 1 of 1 ✔ \[1f/35a1e6\] process \> GTF2BED \[100%\] 1 of 1 ✔ \[fa/b1e933\] process \> TRIM (SRR5655560) \[100%\] 1 of 1 ✔ \[d6/d47ead\] process \> SAM\_FOR\_STRAND (1) \[100%\] 1 of 1, failed: 1 ✘ Error executing process \> 'SAM\_FOR\_STRAND (1)'
>
> Caused by: Process `SAM_FOR_STRAND (1)` terminated with an error exit status (255)
>
> Command executed:
>
> bowtie2 -t -x index.1.bt2 index.2.bt2 index.3.bt2 index.4.bt2 index.rev.1.bt2 index.rev.2.bt2 -1 SRR5655560\_1.atria.fastq -2 SRR5655560\_2.atria.fastq -S "SRR5655560.sam"

The following is simplifies version of the pipeline:

```
// Declare syntax version
nextflow.enable.dsl=2

/*
 * pipeline input parameters
 */

params.reads = "$projectDir/data/*_{1,2}.fastq"
params.ref_genome = "$projectDir/data/Glycine_max.Glycine_max_v2.1.dna.toplevel.fa"
params.gtf_file = "$projectDir/data/Glycine_max.Glycine_max_v2.1.55.gtf"
params.outdir = "results"

log.info """\

        R N A S E Q - N F   P I P E L I N E
        ===================================
        reads            : ${params.reads}
        reference genome : ${params.ref_genome}
        gtf              : ${params.gtf_file}
        outdir           : ${params.outdir}
        """
        .stripIndent()

/*
 * define the `index` process that creates a binary index
 * given the reference genome file
 */

process INDEX {

  input:
  path ref_genome

  output:
  path "index.*"

  script:
  """
  bowtie2-build $ref_genome index
  """
}

/*
 * define the `GTF2BED` process that creates a BED file
 * given the params.gtf_file as input.
 */

process GTF2BED {

    input:
    path gtf

    output:
    path "${gtf.baseName}.bed"

    script:
    """
    gtf2bed --gtf "${gtf}" --bed "${gtf.baseName}.bed"
    """
}

/*
 * define the `TRIM` process that creates trimmed fastq files
 * given the reads fastq file.
 */

process TRIM {

    tag { pair_id }

    input:
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path(trimmed_reads), emit: trimmed_fastqs
    path "${pair_id}/*.atria.log*", emit: logs

    script:
    def (r1, r2) = reads

    trimmed_reads = reads.collect {
        "${pair_id}/${it.baseName}.atria.fastq"
    }

    """
    atria \\
        -r "${r1}" \\
        -R "${r2}" \\
        -o "${pair_id}" \\
        -t 1 \\
        -l 12 \\
        -q 30
    """
}

/*
 * define the `SAM_FOR_STRAND` process that creates a SAM file
 * to be used by `STRAND_CHECK` process to determine
 * the strandedness of each fastq file pair.
 */

process SAM_FOR_STRAND {

    input:
    path index
    tuple val(pair_id), path(reads)

    output:
    path "${pair_id}.sam"

    script:
    """
    bowtie2 -x "${index.baseName}" -1 ${reads[0]} -2 ${reads[1]} -S ${pair_id}.sam
    """
}

workflow {

    index_ch = INDEX(params.ref_genome)

    gtf2bed_ch=GTF2BED(params.gtf_file)

    read_pairs_ch = Channel.fromFilePairs(params.reads, checkIfExists:true)

    trim_out = TRIM(read_pairs_ch)

    sam_for_strand_ch = SAM_FOR_STRAND(index_ch, TRIM.out.trimmed_fastqs)

}

```

I've tried passing the *index* argument directly to the script in the SAM\_FOR\_STRAND Alignmnet process, it gave the "index.1.bt2" does not exist error. Then I tried giving the basename of the index files (as shown in the below). It still does not work, even though the index files are being created properly.

```
 /*
     * define the `SAM_FOR_STRAND` process that creates a SAM file
     * to be used by `STRAND_CHECK` process to determine
     * the strandedness of each fastq file pair.
     */

    process SAM_FOR_STRAND {

        input:
        path index
        tuple val(pair_id), path(reads)

        output:
        path "${pair_id}.sam"

        script:
        """
        bowtie2 -x "${index.baseName}" -1 ${reads[0]} -2 ${reads[1]} -S ${pair_id}.sam
        """
    }

```

Best wishes.

# Answer
When we run Bowtie2, we just need to make sure we supply the index basename using `-x`. This is the index filename prefix minus the `.X.bt2` extension. A safe assumption is that the first file in the collection of output files (sorted in lexicographical order) has the `.1.bt2` suffix. We can then just strip this off to get the required prefix. For example:

```
process SAM_FOR_STRAND {

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id), path(reads)

    output:
    path "${pair_id}.sam"

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 \\
        -x "${idx}" \\
        -1 "${reads[0]}" \\
        -2 "${reads[1]}" \\
        -S "${pair_id}.sam"
    """
}

```
> 1 votes

---
Tags: nextflow
---


# Question
Title: How to handle output from a python script that gives different outputs depending upon the condition - Nextflow
I have created three Nextflow processes: One checks the strandedness of the fastq file and outputs it to a .txt file. Second one executes a python script with different commands depending upon the condition. And the third one aligns them.

```
/*
 * the `CHECK_STRANDEDNESS` process takes SAM file
 * and a BED file
 * and outputs a a txt file containing information about strandedness.
 */

process CHECK_STRANDEDNESS {

    input:
    path "${pair_id}.sam"
    path "${gtf.baseName}.bed"

    output:
    path "${pair_id}_strand_check.txt"

    script:
    """
    infer_experiment.py -i "${pair_id}.sam" -r "${gtf.baseName}.bed" > "${pair_id}_strand_check.txt"
    """
}

/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file
 * and executes commands according to the reads strandedness.
 * The output generated is FASTA files.
 */

process FINDTAIL {

    input:
    path "${pair_id}_strand_check.txt"
    path tuple val(pair_id), path(reads)

    output:

    script:
    """
    #!/usr/bin/env python3

    import subprocess
    import pandas as pd
    import os

    result = pd.read_csv("${pair_id}_strand_check.txt", sep="\r\n", header=None, engine='python')

    failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
    fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
    rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
    fwd_percent = fwd/(fwd+rev)
    rev_percent = rev/(fwd+rev)

    if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
        cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
        cmd2 = "findtail_v1.01 --input_file ${reads[1]}  --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
        print(cmd1)
        subprocess.call(cmd1, shell=True)
        print(cmd2)
        subprocess.call(cmd2, shell=True)

        if fwd_percent > 0.9:
            #Over 90% of reads explained by "1++,1--,2+-,2-+
            #Data is likely FR/fr-secondstrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif rev_percent > 0.9:
            # Over 90% of reads explained by "1+-,1-+,2++,2--
            # Data is likely RF/fr-firststrand
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)

        elif max(fwd_percent, rev_percent) < 0.6:
            #Under 60% of reads explained by one direction
            #Data is likely unstranded
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.sl.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_1.pr.fasta
            cmd3 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.sl.fasta
            cmd4 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > tail_${pair_id}_2.pr.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
            print(cmd3)
            subprocess.call(cmd3, shell=True)
            print(cmd4)
            subprocess.call(cmd4, shell=True)

        else:
            #if strand couldnt be detected
            # we assume it is first-stranded
            # since most illumina reads are
            cmd1 = "findtail_v1.01 --input_file ${reads[0]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_1.fasta
            cmd2 = "findtail_v1.01 --input_file ${reads[1]} --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > tail_${pair_id}_2.fasta
            print(cmd1)
            subprocess.call(cmd1, shell=True)
            print(cmd2)
            subprocess.call(cmd2, shell=True)
        """
}

/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process
 * and try to align them INDIVIDUALLY.
 */

process ALIGN {

    input:
    path index
    path "tail_${pair_id}*.fasta"

    output:
    path "unmapped_tail_${pair_id}*.fasta"

    script:
    """
    bowtie2 -x $index -f "tail_${pair_id}*.fasta" --un "unmapped_tail_${pair_id}*.fasta"
    """
}

```

After creating a *"${pair\_id}\_strand\_check.txt"* file, I input it to the FINDTAIL process that finds reads with polyadenylated ends. It takes fastq files individually, so I have to run it individually for each file.

If fastq file pair is first stranded or second stranded, I executed findtail twice,once for R1 and once for R2, which creates two fasta files. But if it is unstranded, I have to run findtail four times, twice for each file. Thus, depending on the paired-end reads' strandedness, I have different number of inputs. **How do I define the output of this process?**

In the next process I align these files individually to get unmapped reads.

*This issue has been eating my brain away for 5 days. So please let me know if something is unclear.*

Best wishes!

# Answer
You should be able to just use a glob pattern that accounts for all possible output files, regardless of the specified conditions. As long as a single file matches the glob pattern when the process completes (successfully), Nextflow won't complain. I think this should suffice, but please let me know in the comments if it doesn't:

```
output:
tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

```
> 4 votes

---
Tags: nextflow
---


title: umap failed to cluster the cells

text:
# Question
Title: umap failed to cluster the cells
I tried umap visualization with `scanpy` :

```
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata, base=2)
sc.pp.highly_variable_genes(adata, flavor='cell_ranger', n_top_genes=400)

sc.pp.scale(adata, zero_center=True, max_value=None, copy=False, layer=None, obsm=None)
sc.pp.pca(adata, n_comps=50, use_highly_variable=True, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=50)

sc.tl.umap(adata, min_dist=0.5, spread=1.0)
sc.pl.umap(adata, color='fullname', use_raw=False, save='samples_umap.pdf')

```

version

```
anndata     0.7.5
scanpy      1.6.1

```

But the cells can't separate well

I tried another small dataset with `scanpy` using the same parameters as before:

`sc.tl.umap` still failed to down dimension the data properly.

Then I tried the original `umap` package using the same data set:

```
import umap
import umap.plot
mapper = umap.UMAP().fit(adata.X)
umap.plot.points(mapper)

```

Now the original `umap` package can do down dimension very well:

I think there may be something wrong with the `umap` function in `scanpy` Can anyone please let me know the reason? Thanks a lot.

# Answer
> 2 votes
Without more code and data, this can't be debugged. But first question would be does the data pass basic QC (e.g. look at cellranger web summaries or equivalent). If yes, you could trying different normalisation (code not shown) and/or no scaling. Try `sc.tl.umap()` with default params. Plot the PCA, try t-SNE. Plot the basic QC metrics following the introductory scanpy tutorial. That is how I would start debugging

---
Tags: scrnaseq, single-cell, umap
---


title: Read Clustal file in Python

text:
# Question
Title: Read Clustal file in Python
*This question was also asked on StackOverflow*

I have a multiple sequence alignment (MSA) file derived from mafft in clustal format which I want to import into Python and save into a PDF file. I need to import the file and then highlight some specific words. I've tried to simply import the pdf of the MSA but after the highlight command doesn't work.

I need to print the file like this:

`Multi.txt` :

```
CLUSTAL format alignment by MAFFT FFT-NS-i (v7.453)

Consensus       --------------------------------------acgttttcgatatttatgccat
AMP             tttatattttctcctttttatgatggaacaagtctgcgacgttttcgatatttatgccat
                                                      **********************

Consensus       atgtgcatgttgtaaggttgaaagcaaaaatgaggggaaaaaaaatgaggtttttaataa
AMP             atgtgcatgttgtaaggttgaaagcaaaaatgaggggaaaaaaaatgaggtttttaataa
                ************************************************************

Consensus       ctacacatttagaggtctaggaaataaaggagtattaccatggaaatgtatttccctaga
AMP             ctacacatttagaggtctaggaaataaaggagtattaccatggaaatgtaattccctaga
                ************************************************** *********

Consensus       tatgaaatattttcgtgcagttacaacatatgtgaatgaatcaaaatatgaaaaattgaa
AMP             tatgaaatattttcgtgcagttacaacatatgtgaatgaatcaaaatatgaaaaattgaa
                ************************************************************

Consensus       atataagagatgtaaatatttaaacaaagaaactgtggataatgtaaatgatatgcctaa
AMP             atataagagatgtaaatatttaaacaaagaaactgtggataatgtaaatgatatgcctaa
                ************************************************************

Consensus       ttctaaaaaattacaaaatgttgtagttatgggaagaacaaactgggaaagcattccaaa
AMP             ttctaaaaaattacaaaatgttgtagttatgggaagaacaaactgggaaagcattccaaa
                ************************************************************

Consensus       aaaatttaaacctttaagcaataggataaatgttatattgtctagaaccttaaaaaaaga
AMP             aaaatttaaacctttaagcaataggataaatgttatattgtctagaaccttaaaaaaaga
                ************************************************************

Consensus       agattttgatgaagatgtttatatcattaacaaagttgaagatctaatagttttacttgg
AMP             agattttgatgaagatgtttatatcattaacaaagttgaagatctaatagttttacttgg
                ************************************************************

Consensus       gaaattaaattactataaatgttttattataggaggttccgttgtttatcaagaattttt
AMP             gaaattaaattactataaatgttttattataggaggttccgttgtttatcaagaattttt
                ************************************************************

Consensus       agaaaagaaattaataaaaaaaatatattttactagaataaatagtacatatgaatgtga
AMP             agaaaagaaattaataaaaaaaatatattttactagaataaatagtacatatgaatgtga
                ************************************************************

Consensus       tgtattttttccagaaataaatgaaaatgagtatcaaattatttctgttagcgatgtata
AMP             tgtattttttccagaaataaatgaaaatgagtatcaaattatttctgttagcgatgtata
                ************************************************************

Consensus       tactagtaacaatacaacattgga----------------------------------
AMP             tactagtaacaatacaacattggattttatcatttataagaaaacgaataataaaatg
                ************************

```

How can I import the alignment and print in the new PDF with the right alignment of the sequences.

Thanks

# Answer
Ok, figured out a way, not sure its the best one,

nedd to install `fpdf2` (`pip install fpdf2`)

```
from io import StringIO

from Bio import AlignIO  # Biopython 1.80

from fpdf import FPDF # pip install fpdf2

alignment = AlignIO.read("Multi.txt", "clustal")

stri = StringIO()

AlignIO.write(alignment, stri, 'clustal' )

# print(stri.getvalue())

stri_lines = [ i for i in stri.getvalue().split('\n')]

# print(stri_lines)

pdf = FPDF(orientation="P", unit="mm", format="A4")

# Add a page
pdf.add_page()

pdf.add_font('FreeMono', '', 'FreeMono.ttf')

pdf.set_font("FreeMono", size = 8)

for x in stri_lines:
    pdf.cell(0, 5, txt = x, border = 0,  new_x="LMARGIN" , new_y="NEXT",  align = 'L', fill = False)

    # print(len(x))

pdf.output("out.pdf")

```

output pdf `out.pdf` :

Not sure why the file Header is changed, think is something within Biopythion (!!! ???) you can check adding:

```
with open('file_output.txt', 'w') as filehandler:
    AlignIO.write(alignment, filehandler, 'clustal')

```

I had to place a tff FreeMono font (Mono spaced font in my script directory) see: `pdf.add_font('FreeMono', '', 'FreeMono.ttf')` otherwise the alignement won't be printed in the correct way Which fonts have the same width for every character?.

Attached a png of my pdf. See that you can highlight it,

using:

```
pdf.set_fill_color(255, 255, 0)

filling = False
for x in stri_lines:
    if 'Consensus' in x:
        filling = True
    else:
        filling = False
    pdf.cell(0, 5, txt = x, border = 0,  new_x="LMARGIN" , new_y="NEXT",  align = 'L', fill = filling)

```

or something similar you can highlight while printing:
> 2 votes

# Answer
pyMSAviz may help with this problem. pyMSAviz is a Python tool with a convenient CLI/API implementation for visualizing MSA files.

For example, you can output MSA figure in PDF format with the following command.

```
pymsaviz -i Multi.txt -o msa_result.pdf --format clustal \
         --wrap_length 60 --show_count --color_scheme Identity

```

You can also use pyMSAviz's API to configure the appearance of the figure in more detail, such as highlighting.
> 3 votes

# Answer
adding, a second answer because the OP request about

> I need to import the file and then highlight some specific word

was making me uneasy. I kept using the `Pyfpdf2` library and made use of the Biopython `SeqRecord.letter_annotations` attribute to track Nucleotides index composing the words found by RegEx. Tried to use the `page.search_for()` of the `PyMuPDF` library directly on the output pdf of my previous answer but was not getting word accross line breaks \[if anybody know better, please help\].

Here my code, it's a two half problem, first find the word (in my example is '*atat*'), then write everything to pdf, I stole it from from /master/Bio/AlignIO/ClustalIO.py/ def write\_alignment(self, alignment): from class ClustalWriter(SequentialAlignmentWriter) with some modding \[I believe but not absolutely sure that ...A cell by definition is uniform in font and color... \]:

```
import re

from Bio import AlignIO  # Biopython 1.80

from fpdf import FPDF # pip install fpdf2

alignment = AlignIO.read("Multi.txt", "clustal")

word_to_find = 'atat'

for rec in alignment:

    rec.letter_annotations['highlight'] = 'N'*len(rec.seq)

# for rec in alignment:
#     print('rec.id : ', rec.id, 'highlight : ', rec.letter_annotations['highlight'])

for rec in alignment:

    # print(str(rec.seq))

    found = [m.start() for m in re.finditer( '(?='+re.escape(word_to_find)+')' , re.escape(str(rec.seq.replace('-', 'Z'))))]

    # print(rec.id , ' found : ', found)

    for i in found :

            a = list(rec.letter_annotations['highlight'])

            a[i:i+len(word_to_find)] = 'Y'*len(word_to_find)

            rec.letter_annotations['highlight'] = ''.join(i for i in a)

            # print(rec.letter_annotations['highlight'])

max_length = alignment.get_alignment_length()

# stolen from /master/Bio/AlignIO/ClustalIO.py/   def write_alignment(self, alignment): from class ClustalWriter(SequentialAlignmentWriter):
def write_alignment(alignment):

        x = 2

        y = 2

        pdf = FPDF(orientation="P", unit="mm", format="A4")

        # Add a page
        pdf.add_page()

        pdf.add_font('FreeMono', '', 'FreeMono.ttf')

        pdf.set_font("FreeMono", size = 12)

        pdf.set_fill_color(255, 255, 0)

        """Use this to write (another) single alignment to an open file."""
        if len(alignment) == 0:
            raise ValueError("Must have at least one sequence")
        if alignment.get_alignment_length() == 0:
            # This doubles as a check for an alignment object
            raise ValueError("Non-empty sequences are required")

        # Old versions of the parser in Bio.Clustalw used a ._version property
        try:
            version = str(alignment._version)
        except AttributeError:
            version = ""
        if not version:
            version = "1.81"
        if version.startswith("2."):
            # e.g. 2.0.x
            output = f"CLUSTAL {version} multiple sequence alignment"
        else:
            # e.g. 1.81 or 1.83
            output = f"CLUSTAL X ({version}) multiple sequence alignment"

        pdf.cell(0, y, txt = output , border = 0,  new_x='LMARGIN'  , new_y="NEXT" ,  align = 'L', fill = False)
        pdf.cell(0, 4*y, txt = '' , border = 0,  new_x='LMARGIN'  , new_y="NEXT" ,  align = 'L', fill = False)

        cur_char = 0
        # # max_length = len(alignment[0])
        max_length = alignment.get_alignment_length()

        if max_length <= 0:
            raise ValueError("Non-empty sequences are required")

        if "clustal_consensus" in alignment.column_annotations:
            star_info = alignment.column_annotations["clustal_consensus"]
        else:
            try:
                # This was originally stored by Bio.Clustalw as ._star_info
                star_info = alignment._star_info
            except AttributeError:
                star_info = None

        # keep displaying sequences until we reach the end

        rec_cnt = 0

        while cur_char != max_length:
            # calculate the number of sequences to show, which will
            # be less if we are at the end of the sequence
            if (cur_char + 50) > max_length:
                show_num = max_length - cur_char
            else:
                show_num = 50

            # go through all of the records and print out the sequences
            # when we output, we do a nice 80 column output, although this
            # may result in truncation of the ids.
            for record in alignment:

                rec_cnt += 1

                # Make sure we don't get any spaces in the record
                # identifier when output in the file by replacing
                # them with underscores:
                line = record.id[0:30].replace(" ", "_").ljust(50)

                pdf.cell(50, y, txt = line, border = 0,  new_x='RIGHT' , align = 'L', fill = False)

                line_seq = list(str(record.seq[cur_char : (cur_char + show_num)]))

                line_seq_highlight = list(str(record.letter_annotations['highlight'][cur_char : (cur_char + show_num)]))

                for i in range(len(line_seq)):

                    if line_seq_highlight[i] == 'N' :

                        filling = False

                    if line_seq_highlight[i] == 'Y' :

                        filling = True

                    else:

                        filling = False

                    pdf.cell(x, y, txt = line_seq[i], border = 0,  new_x= 'RIGHT' ,  align = 'C', fill = filling)

                pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

                # now we need to print out the star info, if we've got it

                if star_info and rec_cnt == len(alignment):

                    pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

                    rec_cnt = 0

                    # print(star_info)

                    line = (' ' * 50)

                    pdf.cell(50, y, txt = line, border = 0,  new_x='RIGHT' , align = 'L', fill = False)

                    star_info_seq = list(str(star_info[cur_char : (cur_char + show_num)]))

                    for i in range(len(star_info_seq)):

                        pdf.cell(x, y, txt = star_info_seq[i], border = 0,  new_x= 'RIGHT' ,  align = 'C', fill = False)

                pdf.cell(x, 2*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

            cur_char += show_num

            pdf.cell(x, 3*y, txt = ' ' , border = 0,  new_x="LMARGIN" , new_y="NEXT" ,  align = 'L', fill = False)

        pdf.output("out.pdf")

write_alignment(alignment)

```

I tested it, with the `Multi.txt` Clustal alignment from above, output is

`out.pdf` :

.

One last note:

needed to replace `-` with `Z` in the Consensus string to be able to get reliable results from the RegEX line : `found = [m.start() for m in re.finditer( '(?='+re.escape(word_to_find)+')' , re.escape(str(rec.seq.replace('-', 'Z'))))]` , not sure what is going on here, got the snippet from a SO post.
> 2 votes

---
Tags: python, sequence-alignment, phylogenetics, multiple-sequence-alignment
---


# Question
Title: Reference for the statistical methodology used in joint genotype calling
I am searching for a good reference that explains the statistical methodology used in what is variously known as Joint Genotyping, Cohort Genotyping or Multi-sample Genotyping. One with a pseudo code implementation would also be ideal.

I have found several discussing the approach in very general terms and the relative benefits of the technique but nothing with actual equations or algorithms. Failing that, if there is an open source implementation that people would consider to be a reference implementation and that has excellent code documentation that would also be great.

# Answer
> 2 votes
To the best of my knowledge, the first joint calling model was published by Le & Durbin (2011) and refined in Li (2012). They come with detailed derivation. Both samtools and older GATK are using this model. GATK now uses a different model which only works with high-coverage data but is faster to compute. There are other alternatives (e.g. in freebayes).

# Answer
> 2 votes
Here is a recent example of the derivation of the genotype model from Octopus genotype calling paper:

> All samples have known ploidy and copy number, so the likelihood function of reads $R$ given genotype $g$ is:
>
> $p(R|g) = \mathop {\prod }\limits\_{n = 1}^{|R|} \frac{1}{{|g|}}\mathop {\sum }\limits\_{i = 1}^{|g|} p(r\_n|g\_i)$
>
> where $|g|$ is the ploidy and $|R|$ is the number of reads. The joint genotype posterior for $S$ samples is therefore given by:
>
> $p({\bf{g}}|{{R}},{\cal{M}}\_{{g}}) \propto p({\bf{g}}|{\cal{M}}\_g)\mathop {\prod }\limits\_{s = 1}^S p({{R}}\_s|{\bf{g}}\_s)$
>
> where the genotype prior model, $\cal{M}$, is either the uniform or HWE-coalescent prior. Unfortunately, the number of genotype combinations g grows exponentially in the number of samples S, so we cannot evaluate the full posterior distribution in general. Therefore, other than for trivial cases, we first approximate the sample marginal genotype posterior distribution under the HWE model (without mutations) $(|,\cal{M})$ and use these marginal probabilities to select K genotype combinations $1,…,$ (K is user-defined) to evaluate under the full joint genotype model. The approximate posterior marginals are computed with expectation maximization.

The source code can be found here.

---
Tags: genotyping
---


# Question
Title: How to manually trim FASTA file sequences with the information provided in the header and store it into a new FASTA file? - Python
I have a FASTA file that has reads like this:

> \>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) CttttatttgttGTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

As you can see, there is information about the start "`start=2`" and end "`end=12`" of a sequence within the header. I would like to slice the sequence like \[start:end\] and keep the rest of it. e.g

The part I would like to trim:

> Ctttggtttcctttt

And after trimming, I would like to keep the rest of the read:

> GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

So the new FASTA file will have reads like this:

> \>SRR5655563.745 745 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%) GTGTGTGAGGTTTGATTTGATGGGAAAATATCTTGAATCTGCGGCGAGGTTGGAAGAGCTATCGCGGATTGTGTCATCTGCTGCGAAGCCCAATAGGTCAAAGGGAACGCTACC

I have written a python script that takes a FASTA file as an input and outputs a trimmed FASTA file, like this:

```
tailtrim.py --fasta input.fasta --output output.tail.trimmed.fasta

```

The code does trim the reads where I want it to.

```
import numpy as np
import argparse
import re

def main():
    parser = argparse.ArgumentParser()

    parser.add_argument("--fasta", type=str, required=True, help="input fasta file")
    parser.add_argument("--output", type=str, required=True, help="output tail.trimmed.fasta file")

    args = parser.parse_args()

    out_fasta = open(args.tail.trimmed.fasta, 'w')

    for line in open(args.fasta):
        header = line.rstrip()
        read = next(args.fasta).rstrip()

    #Extract header into a list
        header_lst = []
        header_lst= header.split(" ")

        tmp_start_lst = []
        tmp_end_lst = []
        tmp_lst = []

        for i in header_lst:
            if i.startswith("(") and i.endswith(")"):
                tmp_lst.append(i.split(","))
        #print(tmp_lst)

    # Extract the start and end numbers in a list
        for i in tmp_lst:
            for str in i:
                if str.startswith("start="):
                    tmp_start_lst.append(str)
                if str.startswith("end="):
                    tmp_end_lst.append(str)

        pattern = re.compile(r'\d+$')

        start_lst = []
        for str in tmp_start_lst:
            match = re.search(pattern, str)
            start_lst.append(int(match.group()))

        end_lst = []
        for str in tmp_end_lst:
            match = re.search(pattern, str)
            end_lst.append(int(match.group()))

        #Extract reads into a list to manipulate i.e. trim
        read_lst = []
        read_lst.append(read)

        new_lst = []

        for str in read_lst:
            for num_start in start_lst:
                for num_end in end_lst:
                    s = str.replace(str[(num_start-1):num_end],"")
                    new_lst.append(s)

```

I do get the trimmed reads, but what I'm having trouble figuring out is how to write these trimmed reads with their corresponding headers to a new FASTA file.

Also if there is any way that I could improve my code or employ better strategy to approach the problem, I would love to know.

Best wishes.

Edit: I tried running the solution by the super helpful user @Steve.

The command I ran:

```
python3 tailTrimming.py --fasta tailtrim.fasta --output testtrim.trim.fasta

```

The command line error:

```
Traceback (most recent call last):
  File "/home/x/Desktop/fastq/tailTrimming.py", line 61, in <module>
    main()
  File "/home/x/Desktop/fastq/tailTrimming.py", line 58, in main
    SeqIO.write(trim_reads(fasta), out, "fasta")
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/__init__.py", line 478, in write
    count = writer_class(fp).write_file(sequences)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 209, in write_file
    count = self.write_records(records)
  File "/home/x/.local/lib/python3.10/site-packages/Bio/SeqIO/Interfaces.py", line 193, in write_records
    for record in records:
  File "/home/x/Desktop/fastq/tailTrimming.py", line 20, in trim_reads
    result = match.groupdict()
AttributeError: 'NoneType' object has no attribute 'groupdict'

```

This is the whole input file I supplied:

```
>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
tttctttctttAGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATATTGTATTGTACCCTCAGATTTCTCGCACAAT
>SRR5655563.10087 10087 length=126 (type=T,start=2,end=16,length=15,identity=80%)
GtttgtttcattttttCCAAGGAAGAAATTCGACCTTCCAGCGATTTGTGGAGGTCATTCATCTGTTGAGTAAAGGATTGGTGCGCGACCTTGTAATCCTTGTGGAGGGTGTAG
>SRR5655563.13902 13902 length=126 (type=T,start=2,end=16,length=15,identity=80%)
CtttttcttcttgtttGATGTAAAAAATCTTCGGTATCTAAAACGTATGAAACCTTAGATGAACTGCAAGCTATCTTACACTCTAGACCATTGATAATTT
>SRR5655563.18489 18489 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttttcttcttCTTGGAAGTTTTGGATTCTGCAGGGACATTTTCTTTCTTTTCTCCATCTCCATCAGTCTCAACACCTGTCTTATCCTGAGCACCTTGAGACTCATCTTGACCATCTTTACTTT
>SRR5655563.20782 20782 length=126 (type=T,start=2,end=15,length=14,identity=71.4286%)
CttcattttagttttCAACTCCTCTACCAACTGGGATGGTGCATTCAACTGTAGAACTGTTTCACGAATCTCCTTGAGAACACTGAAGTCTCCCTCAGTTAACTTCTTTTTCAAATTATTGA
>SRR5655563.23972 23972 length=126 (type=T,start=2,end=29,length=28,identity=75%)
CttttttgttgttctcttctttcttctttGCATGGCTTGATGATGACTGATCCGATTGATTATTTCTTGATAGAATATCAATGTCATCACCAATTGTTGAGT
>SRR5655563.31739 31739 length=126 (type=T,start=1,end=10,length=10,identity=80%)
ttttcttcttCTTGCATGAAAGTGGGTTTGGACCCTTAGCTCTTTTTCTCTTGAACTGAGGTTT
>SRR5655563.34244 34244 length=126 (type=T,start=2,end=14,length=13,identity=76.9231%)
CtttcctttattttCTTGGAGGTTAGCCAGGATATTGGGATCAAAGCATGTAGCAAAACACACCTTGCTATTTCTTGCTCGTACAGATACATGGGATAGGA
>SRR5655563.38334 38334 length=126 (type=T,start=2,end=12,length=11,identity=81.8182%)
CtttattgttttCTCATCACCCATATTATATTTCTTAAGGGCAGATTCCAACTCTTCCATTGTGATGTACCCGCTCTTATCATTGTCAAAATATTCAAAGGCTTTATATAGGCGATCCTCTCTTTC

```

And this is the whole output I have gotten:

```
>SRR5655563.9676 9676 length=126 (type=T,start=1,end=11,length=11,identity=81.8182%)
AGTTCCAATAGTCTTGACCACCTTTTATCAACAGAAGCTTGGCCTTGTCCAGTTATATAT
TGTATTGTACCCTCAGATTTCTCGCACAAT

```

Just one read, that's all. Whereas I would like to keep all the sequences because all of them have this `start`, `end` attribute. The sequence being on two lines in the output is not something I did. It is like this in the output file. So it's no longer FASTA file, right?

Thanks in advance for any help!

# Answer
Here's one way using Biopython and the SeqIO interface to read and write SeqRecord objects. Using a generator, we can produce (yield) trimmed sequence records that can then be written to the output handle. This avoids storing an entire list of trimmed sequences in memory and is usually the way to go, assuming you're dealing with large files containing thousands of records:

```
import argparse
import re

from Bio import SeqIO

```

```
def trim_reads(fileobj):

    pattern = (
        r'type=(?P<type>\w+),'
        r'start=(?P<start>\d+),'
        r'end=(?P<end>\d+),'
        r'length=(?P<length>\d+),'
        r'identity=(?P<identity>\d+(?:\.\d+)?)%'
    )

    for seq_record in SeqIO.parse(fileobj, "fasta"):
        match = re.search(pattern, seq_record.description)
        result = match.groupdict()

        start = int(result['start']) - 1
        end = int(result['end'])

        head = seq_record.seq[:start]
        tail = seq_record.seq[end:]

        seq_record.seq = head + tail

        yield seq_record

```

```
def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "--fasta",
         type=str,
         required=True,
         help="input fasta file",
    )
    parser.add_argument(
        "--output",
        type=str,
        required=True,
        help="output tail.trimmed.fasta file",
    )

    return parser

```

```
def main():

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

    with open(args.fasta) as fasta, open(args.output, "w") as out:
        SeqIO.write(trim_reads(fasta), out, "fasta-2line")

```

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

```
> 3 votes

---
Tags: python, fasta
---


# Question
Title: Mismatch between mutant gene and reference gene outside the site of mutation
I am tasked with designing primers for a particular mutant (Target Gene Locus: At1g28490 i.e. SYP61) of Arabidopsis, obtained from Gabi-Kat:

* The NCBI reference sequence for the genome goes: GGCTTTATAAGATTCTCCAGTAGCCCAAAAATAAACTTTCAGCGAATCTCAAAGCCAACCGCAGGTGATGCGTTTTGCTATGTAAGAGAGACGAATCTCACGCCGGAATCTCTTTTACTTCACCGTCGTTTGTTCTCCGTCGCGTCACCTTTCGATTAATTCCTTCACAGGGTAGCAAAATTTTGTTGATTTCGGG$\color{red}{\textrm{\[}}$ATGTCTTCAGCTCAAGATCCATTCTACATTGTTAAAGAAGAGATTCAGGATTCTATTGATAAGTTGCAATCTACATTCCACAAATGGG**AGCG**T**ATCTCTCCCGACATGG**G**TGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGT**A͟G͟A͟T͟G͟A͟G͟C͟T͟G͟G͟A͟A͟A͟A͟A͟G͟C͟G͟A͟T͟T͟A͟C͟T͟G͟T͟T͟G͟C͟A͟G͟C͟A͟A͟A͟A͟G͟A͟T͟C͟C͟T͟T͟C͟ATGGTATGGCATTGATGAAGCTGAGCTTGAAAAACGGAGGAGATGGACTAGTAATGCTAGGACACAGGTGCGGAATGTGAAGTCTGGTGTTTTAGCTGGAAAGGTTAGTTCTGGAGCAGGTCATGCAAGTGAAGTGCGTAGAGAACTAATGAGAATGCCAAACTCGGGTGAAGCAAGTAGATATGATCAGTATGGAGGGAGAGACGATGATGGTTTTGTACAGTCAGAATCAGATAGGCAAATGCTGTTGATAAAGCAACAAGACGAAGAATTGGATGAGCTTAGTAAAAGTGTACAGAGAATTGGAGGAGTGGGGCTCACTATACACGACGAACTCGTTGCACAGGAGAGAATAATAGATGAATTGGATACGGAGATGGACAGCACAAAGAACCGTCTTGAATTTGTTCAGAAAAAAGTGGGAATGGTAATGAAGAAAGCTGGAGCAAAAGGTCAAATGATGATGATATGTTTCTTGCTCGTCTTGTTCATCATCCTATTCGTTCTTGTCTTCTTGACCTAA$\color{red}{\textrm{\]}}$AGAACAGAGCCGCGGGATCTTATAGATGGGAGTTTAATGTCTTACATTATCATGTTTGTGGTGAAACAAATGGAATATTTACATAAGGTTGTCTGTTTTTTCTTCTCCTATCCAATCCCACATATAATTTTGAACTTTGGCTGTGTGATCTCTAGAGAACCTGAAATTAAATACTTTCATAGCTTGATCCCATGGCTAAGTTCCT
* Gabi-Kat provides the following sequence from the mutant as a confirmatory match with the genome locus: **AGCG**G**ATCTCTCCCGACATGG**C**TGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGT**C͟C͟T͟C͟C͟T͟A͟T͟T͟C͟A͟C͟C͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟C͟C͟C͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟

I have **bold-faced** the common region to both sequences. What explains the *complete mismatch* of the trailing part (u͟n͟d͟e͟r͟l͟i͟n͟e͟d͟) of the Gabi-Kat sequence with the ref genome?

It is **not** the case that the mutation happened there; we are in the coding region (demarcated by red square brackets in the NCBI sequence) while Gabi-Kat's mutation is in an untranslated region.

## Update 1:

Running Ensembl Blast on the GabiKat seq returns a 92.9% identicality with SYP61. Their sequence for the gene (taken from TAIR10.1) goes:

AGCGGATCTCTCCCGACATGGCTGATCAAGCACATGTGGCGAAGGAGCTTGTTGCTACTTGTGGAAGCATTGAGTGGCAGGTA͟C͟T͟C͟G͟T͟A͟T͟T͟T͟A͟G͟T͟C͟A͟C͟G͟G͟T͟T͟C͟T͟A͟T͟T͟A͟A͟T͟C͟T͟C͟G͟T͟C͟T͟T͟T͟A͟T͟A͟T͟A͟T͟A͟A͟

Note the difference in the trailing part from NCBI ref; this almost coincides with the Gabi-Kat seq!

So, is the NCBI ref-genome and the data for the particular TAIR genome inaccurate? Maybe, TAIR updated their whole genome set but forgot to update data for individual genes?

## Update 2:

My proposed explanation for the anomaly makes no sense. TAIR 7 (2004) has the same sequence. It looks like I am missing out on something more fundamental.

# Answer
> 0 votes
**Conclusion** The evidence suggests the region you are looking at **Gabi-kat** appears either 27045 downstream or 192250 upstream of syntaxin SNARE (SYP61), but is definitely not within the gene. It is chance homology (e.g. a pseudo-gene) that appears to occur within syntaxin (SYP61). The intergenic region could be much closer to the SYP61 gene if the same reference genome was used for both analyses Genbank LR782542.1 versus Genbank NC\_003070.9 (different reference genomes in my analysis). Hence, the University is probably correct.

What I suggest is the genome of *A. thaliana* is downloaded and Blast analysis is performed locally. This way differences in the reference genome are removed.

*The analysis is as follows*:

The target protein in question is involved in vesicle trafficking (here) called SNARE, comprising coiled-coil domains and are highly conserved, for example show 100% conservation between *Arabidopsis thaliana*, *Arabidopsis suecica* and *Camelina sativa*

> MSSAQDPFYIVKEEIQDSIDKLQSTFHKWERISPDMGDQAHVAKELVATCGSIEWQVDELEKAITVAAKDP

If Gabi-kat was a protein its sequence would be ...

> RISPDMADQAHVAKELVATCGSIEWQVLLFTHGSINPLFIY

Here's a blast of the whole regions of the Gabi-kat. The region in green is the SNARE proteins ...

So the issue is what is the terminal 3' protein and Blast that, i.e. the final bit that matches no known plant and the locus you flagged as being anomalous ..

This is basically zero result via Blast, the protein with low homology are lower metazoans (animals) and a few bacteria ... a long way from plants. In other words this isn't a known protein of natural origin at least on NCBI.

The locus is literally like no protein on earth.

However ... Blasting the nucleotide ...

1. Reference *Arabidopsis thaliana* = syntaxin (SNARE; SYP61)
2. Gabi-Kat *Arabidopsis thaliana* top two Blast hits are ...

* Chromosome 1 at 10039082 - 10039208, or 9819787 to 981991 or
* Both value 1e-42
* Both dentity is 93%

So Gabi-Kat is a real part of the *Arabidopsis thaliana*, **Gabi-kat is a non-coding region**. *Thats the major conclusion*.

The question where is the location of SYP61? And that is

* Chromosome 1 at 10012037 to 10020311

Described here from NCBI's full report

Thus the mutant is 27045bp downstream or 192250 bp upstream **IF** the same reference genome is used for the blasting (Genome LR782542.1) and the Genbank report NC\_003070.9. They could easily be much closer to each other is the same reference genome was used for both sets of analysis.

---

To answer the comments

* I \[the OP\] do not follow your \[my\] blast analysis
* My \[the OP\] blasts gives ... \['hits within the gene'\]
* etc ....

It appears to me there is a psuedo-gene phenomenon occurring which is making you think a mutant is located within a gene. However, if you disagree thats perfectly okay and the original version of this answer I had assumed the Gaba-kat team had got this wrong (I always from the basis of supporting the OP).

NCBI was the database repository in the remote analysis, but there's an issue with multiple references ... so to reiterate I suggest you perform this analysis as a local blast to remove multiple reference phenomena.

To conclude I personally don't think there is a basis to believe the University responsible for Gaba-kat has got this wrong https://www.gabi-kat.de. However, if you feel otherwise then I would recommend contacting them with your analysis and any supportive analysis for example from this thread or elsewhere. Please note, I don't have a connection with this university.

---
Tags: gene, genomics, blastn, blastp
---


# Question
Title: How to run a nextflow process for each file generated by another process separately - Input tuple does not match input set cardinality error
I have a process called `FINDTAIL` that generates different number of files depending on the input data. Its either 2 files or four files i.e.

> `1. read_1{,.pr,.sl}.fasta and read_2{,.pr,.sl}.fasta`

or

> `2. read_1.pr.fasta, read_1.sl.fasta, read_2.pr.fasta and read_2.sl.fasta`

```
/*
 * the `FINDTAIL` process takes "${pair_id}_strand_check.txt" file
 * and executes commands according to the reads strandedness.
 * The output generated is FASTA files.
 */

process FINDTAIL {

    tag { pair_id }

    debug true

    input:
    path strand_check_file
    tuple val(pair_id), path(reads)

    output:
    tuple val(pair_id), path("tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

    script:
    """
#!/usr/bin/env python3

import subprocess
import pandas as pd
import os

result = pd.read_csv("${strand_check_file}", sep="\\r\\n", header=None, engine='python')

failed = float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', ''))
fwd = float(result.iloc[2,0].replace('Fraction of reads explained by "1++,1--,2+-,2-+": ', ''))
rev = float(result.iloc[3,0].replace('Fraction of reads explained by "1+-,1-+,2++,2--": ', ''))
fwd_percent = fwd/(fwd+rev)
rev_percent = rev/(fwd+rev)

if float(result.iloc[1,0].replace('Fraction of reads failed to determine: ', '')) > 0.50:
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

if fwd_percent > 0.9:
    #Over 90% of reads explained by "1++,1--,2+-,2-+
    #Data is likely FR/fr-secondstrand
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_1.pr.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_2.pr.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

elif rev_percent > 0.9:
    # Over 90% of reads explained by "1+-,1-+,2++,2--
    # Data is likely RF/fr-firststrand
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)

elif max(fwd_percent, rev_percent) < 0.6:
    #Under 60% of reads explained by one direction
    #Data is likely unstranded
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_1.pr.fasta"
    cmd3 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    cmd4 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type pr --d > " + "tail_${pair_id}_2.pr.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)
    print(cmd3)
    subprocess.call(cmd3, shell=True)
    print(cmd4)
    subprocess.call(cmd4, shell=True)

else:
    #if strand couldnt be detected
    # we assume it is first-stranded
    # since most illumina reads are
    cmd1 = "findtail_v1.01 --input_file " + "${reads[0]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_1.sl.fasta"
    cmd2 = "findtail_v1.01 --input_file " + "${reads[1]}" + " --seqlength 500 --endgap 2 --taillength 10 --identity 10 --ptype A --stype T --output_format fasta --output_type sl --d > " + "tail_${pair_id}_2.sl.fasta"
    print(cmd1)
    subprocess.call(cmd1, shell=True)
    print(cmd2)
    subprocess.call(cmd2, shell=True)
    """
}

```

I want to pass these `fasta file` outputs of the `FINDTAIL` process into an `ALIGN` process that runs on each file separately.

The `ALIGN` process:

```
/*
 * define the `ALIGN` process that outputs unmapped files in FASTA format
 * It will take the FASTA files creates by the `FINDTAIL` process
 * and try to align them INDIVIDUALLY.
 */

fasta_file_ch = find_tail_ch.flatten()

process ALIGN {

    tag { pair_id }

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id)
    file(fasta_file) from fasta_file_ch

    output:
    tuple val(pair_id), path("unmapped_tail_${pair_id}_{1,2}{,.pr,.sl}.fasta")

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 -x "${idx}" -f "${fasta_file}" --un "unmapped_tail_${pair_id}/${it.baseName}.fasta"
    """
}

```

I have tried using the `flatten()` function to format my `find_tail_ch` but it still does not work.

Best wishes~

# Answer
> 4 votes
If the output channel produces a tuple where the first element is a simple value and the second element is a collection of files, like in this case, you can simply use the transpose operator to produce a transposition where the first element in each resulting tuple is the *pair\_id* and the second is a FASTA file. You'll then need to update your process definition to receive this tuple, for example:

```
process ALIGN {

    tag { pair_id }

    input:
    path index_files, stageAs: "bt2_index/*"
    tuple val(pair_id), path(fasta_file)

    output:
    tuple val(pair_id), path("${fasta_file.baseName}.unmapped.fasta")

    script:
    def idx = index_files[0].getBaseName(2)

    """
    bowtie2 \\
        -x "${idx}" \\
        -f "${fasta_file}" \\
        --un "${fasta_file.baseName}.unmapped.fasta"
    """
}

```

Your workflow block might look like:

```
workflow {

    bt2_index = ...

    ...

    align_input_ch = find_tail_ch.transpose()

    ALIGN( bt2_index, align_input_ch )
}

```

---
Tags: nextflow
---


# Question
Title: Availability of information on genes in Gnomad VCF data
I was wondering does the gnomad genome data that is downloaded in the vcf format on variants contain information on what is the nearest gene and is the genomic location available exome/intron?

if yes then what are the columns that contain this data? Is it the INFO column?

I have a list of SNPs with chromosome position and the reference and alternative allele information and for each, I would like to know what is he nearest gene and the location:exome/intron.

example data

# Answer
> 4 votes
This question has gone through some iterations and commentary, so I've decided to just answer it, though it's more of a "homework" question. What you need to do is to look at the VCF file, the main thing you want is already in there.

Go to the VCF of interest. I chose the mtDNA because it's smaller and more convenient. You can `zless` into the bgz file to see what it looks like, but if you want uncompressed you can run this (head -1000 just to be a tractable example):

```
gzip -dc gnomad.genomes.v3.1.sites.chrM.vcf.bgz | grep -v intergenic | head -1000 > mt.vcf

```

I `grep -v intergenic` because one of the first things that pops up is that most of the first variants in the file are specifically annotated as `intergenic_variant`, so that answers part of your question. For those variants you can just run a bedtools closest against the right GFF3 version, if you think the nearest gene will be informative.

For variants in genes, what you want is already in the VCF. See the entry for this variant:

```
chrM    1492    .       A       C       .       npg     filters=npg;variant_colla
psed=A1492C;vep=C|non_coding_transcript_exon_variant|MODIFIER|MT-RNR1|ENSG0000021
1459|Transcript|ENST00000389680|Mt_rRNA|1/1||ENST00000389680.2:n.845A>C||845|||||
1||1|SNV||HGNC|HGNC:7470|YES||||||||1||||||||||||;base_qual_hist=0|0|0|0|0|0|0|0|
0|0;position_hist=0|0|0|0|0|0|0|0|0|0;strand_bias_hist=0|0|0|0|0|0|0|0|0|0;weak_e
vidence_hist=0|0|0|0|0|0|0|0|0|0;contamination_hist=0|0|0|0|0|0|0|0|0|0;heteropla
smy_below_min_het_threshold_hist=2|0|0|0|0|0|0|0|0|0;excluded_AC=2;AN=56432;AC_ho
m=0;AC_het=0;hl_hist=0|0|0|0|0|0|0|0|0|0;dp_mean=2976.53;mq_mean=.;tlod_mean=.;AF
_hom=0.00000;AF_het=0.00000;max_hl=0.00000;hap_AN=2680|1537|868|603|34|282|91|147
83|701|934|3144|2732|663|2977|4724|5672|126|1|1298|366|7|393|3079|6037|1234|819|5
46|12|89;hap_AC_het=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0;hap
_AC_hom=0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0|0;hap_AF_hom=0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00|0.00000e+00;hap_AF_het=0.00000e+00|0.00000e+00|0.00 ...
# this is of course a very long line!!! it goes on but all the interesting stuff already happened.

```

You can see that this variant lies within the *RNR1* rRNA gene of the mitochondrial genome (`MT-RNR1` in the file). Granted it mangles it a little bit by describing a variant in a non-protein-coding gene as a `non_coding_transcript_exon_variant`, but that's the kind of thing that you can work around.

Indeed, if you look at the header of the file, you can see that quite extensive annotations are included (where available):

```
##INFO=<ID=vep,Number=.,Type=String,Description="Consequence annotations from Ensembl VEP; note that the SINGLE_EXON flag and END_TRUNC filters have been removed from the LOFTEE annotations to avoid misinterpretation in context of the mitochondrial genome. Format: Allele|Consequence|IMPACT|SYMBOL|Gene|Feature_type|Feature|BIOTYPE|EXON|INTRON|HGVSc|HGVSp|cDNA_position|CDS_position|Protein_position|Amino_acids|Codons|ALLELE_NUM|DISTANCE|STRAND|VARIANT_CLASS|MINIMISED|SYMBOL_SOURCE|HGNC_ID|CANONICAL|TSL|APPRIS|CCDS|ENSP|SWISSPROT|TREMBL|UNIPARC|GENE_PHENO|SIFT|PolyPhen|DOMAINS|HGVS_OFFSET|MOTIF_NAME|MOTIF_POS|HIGH_INF_POS|MOTIF_SCORE_CHANGE|LoF|LoF_filter|LoF_flags|LoF_info">

```

In conclusion, **look at your data**!

---
Tags: vcf, genome, gene, exon, gnomad
---


title: VCF: false insertions?

text:
# Question
Title: VCF: false insertions?
I'm working with gvcf files, containing also non-variant positions. I have noticed some insertions that I'm not able to correctly interpret or handle. An example (selected samples from a huge vcf):

```
NC_013991.2 1403    .   G   GAA 176,32  .   AC=4;AF=0.014;AN=286;DP=1294;ExcessHet=0.0001;FS=0;InbreedingCoeff=0.5011;MLEAC=3;MLEAF=0.01;MQ=28.09;QD=28.2;SOR=3.912 GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:1403_G_GAA:90,6,0:1403  0/0:2,0:2:6:.:.:0,6,88:.    ./.:0,0:1:.:.:.:0,0,0:. ./.:0,0:0:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 0/0:2,0:2:3:.:.:0,3,45:.    ./.:1,0:1:.:.:.:0,0,0:. ./.:1,0:1:.:.:.:0,0,0:. 1/1:0,3:3:9:.:.:123,9,0:.   0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:2,0:2:6:.:.:0,6,90:.    0/0:9,0:9:21:.:.:0,21,315:.
NC_013991.2 1404    .   A   .   .   .   DP=1286 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21
NC_013991.2 1405    .   A   .   .   .   DP=1282 GT:AD:DP:RGQ    0/0:2:2:6   0/0:2:2:6   0/0:2:2:6   ./.:0:0:0   0/0:1:1:3   0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:3:3:9   0/0:1:1:3   0/0:1:1:3   0/0:2:2:6   0/0:9:9:21

```

Here, the reference sequence from position 1403 to 1405 is "GAA". At site 1403, an alternative allele "GAA" is indicated, which I would interpret as an insertion of "AA". Two individuals thus would have genotype "GAAAA" here while the rest has "GAA".

I've seen multiple cases like this, here with longer alternative allele (again, selected samples, not necessarily the same as before):

```
NC_013991.2 23  .   T   TAATACTTATGTGTTG    66,98   .   AC=2;AF=0.009346;AN=214;DP=314;ExcessHet=0.0105;FS=0;InbreedingCoeff=0.3797;MLEAC=3;MLEAF=0.014;MQ=51.77;QD=33.49;SOR=2.303 GT:AD:DP:GQ:PL  1/1:0,2:2:6:85,6,0  0/0:3,0:3:9:0,9,135 0/0:3,0:3:9:0,9,95  0/0:1,0:1:3:0,3,45  ./.:1,1:2:.:0,0,0   0/0:3,1:4:0:0,0,102
NC_013991.2 24  .   A   .   .   .   DP=322  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:2:2:6   0/0:4:4:12
NC_013991.2 25  .   A   .   .   .   DP=339  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:4:4:12
NC_013991.2 26  .   T   .   .   .   DP=353  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 27  .   A   .   .   .   DP=361  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 28  .   C   .   .   .   DP=368  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:15
NC_013991.2 29  .   T   .   .   .   DP=386  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 30  .   T   .   .   .   DP=388  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 31  .   A   .   .   .   DP=398  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 32  .   T   .   .   .   DP=405  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 33  .   G   .   .   .   DP=419  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   ./.:4:5:0
NC_013991.2 34  .   T   .   .   .   DP=427  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 35  .   G   .   0,09    LowQual BaseQRankSum=-1;DP=433;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=51.42;MQRankSum=1;ReadPosRankSum=0    GT:DP:RGQ   0/0:3:9 0/0:3:9 0/0:3:9 0/0:1:3 0/0:3:9 0/0:5:0
NC_013991.2 36  .   T   .   .   .   DP=441  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:1:1:3   0/0:3:3:9   0/0:5:5:12
NC_013991.2 37  .   T   .   .   .   DP=447  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   0/0:5:5:12
NC_013991.2 38  .   G   .   .   .   DP=449  GT:AD:DP:RGQ    0/0:3:3:9   0/0:3:3:9   0/0:3:3:9   0/0:2:2:6   0/0:3:3:9   ./.:4:5:0

```

The alternative allele "TAATACTTATGTGTTG" is exactly the reference. Of course, the sample in question is homozygous. It seems to me these are "false" insertions: they look like insertions, but aren't. I can write a simple awk script trimming the the alternative allele, but there are cases where there are real insertions, and I don't know how to distinguish those.

I've tried the solution mentioned in this answer (double pass through bcftools norm) but it doesn't resolve these issues.

I could spend a few days coding a C++ program to do this, but surely someone with better understanding of vcf, C++, etc, has already done this?

Here's a case I really don't know how to interpret:

```
NC_013991.2 3321    .   A   ATTTCATT    638,04  .   AC=14;AF=0.044;AN=320;DP=1873;ExcessHet=0;FS=0;InbreedingCoeff=0.5521;MLEAC=10;MLEAF=0.031;MQ=29.5;QD=34.42;SOR=1.739   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:2,0:2:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,16:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:33,0:33:96:.:.:0,96,1440:.  0/0:15,0:15:39:.:.:0,39,585:.
NC_013991.2 3322    .   T   .   .   .   DP=1870 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:2:2:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:33:33:96    0/0:14:14:36
NC_013991.2 3323    .   T   .   .   .   DP=1859 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   ./.:1:2:0   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:36
NC_013991.2 3324    .   T   .   .   .   DP=1849 GT:AD:DP:RGQ    0/0:2:2:6   0/0:4:4:12  0/0:2:2:6   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:1:1:3   0/0:8:8:21  0/0:34:34:99    0/0:14:14:33
NC_013991.2 3325    .   TGG T   639,64  .   AC=14;AF=0.044;AN=318;DP=1843;ExcessHet=0;FS=0;InbreedingCoeff=0.5416;MLEAC=11;MLEAF=0.035;MQ=29.53;QD=28.53;SOR=1.739  GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:21:.:.:0,21,315:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3326    .   G   T   1294,06 .   AC=15;AF=0.055;AN=274;BaseQRankSum=0;DP=1833;ExcessHet=0;FS=2.451;InbreedingCoeff=0.5364;MLEAC=17;MLEAF=0.062;MQ=56.66;MQRankSum=0;QD=16.38;ReadPosRankSum=0;SOR=0.828  GT:AD:DP:GQ:PGT:PID:PL:PS   ./.:0,2:2:.:.:.:0,0,0:. ./.:0,4:4:.:.:.:0,0,0:. ./.:0,2:2:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 1|1:0,1:1:3:1|1:3309_G_A:45,3,0:3309    1/1:0,4:4:12:.:.:145,12,0:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3327    .   G   .   0,02    LowQual DP=1832;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=29   GT:DP:RGQ   0/0:2:0 0/0:4:0 0/0:2:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:0 0/0:1:3 0/0:8:18    0/0:34:99   0/0:14:33
NC_013991.2 3328    .   G   A   650,66  .   AC=15;AF=0.048;AN=312;BaseQRankSum=0;DP=1826;ExcessHet=0;FS=3.31;InbreedingCoeff=0.5371;MLEAC=12;MLEAF=0.038;MQ=34.34;MQRankSum=0;QD=27.23;ReadPosRankSum=0;SOR=0.324   GT:AD:DP:GQ:PGT:PID:PL:PS   1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 1|1:0,4:4:12:1|1:3321_A_ATTTCATT:180,12,0:3321  1|1:0,2:2:6:1|1:3321_A_ATTTCATT:90,6,0:3321 ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. ./.:0,1:1:.:.:.:0,0,0:. 0/0:1,0:1:3:.:.:0,3,45:.    0/0:1,0:1:3:.:.:0,3,25:.    0/0:8,0:8:18:.:.:0,18,270:. 0/0:34,0:34:99:.:.:0,99,1485:.  0/0:14,0:14:33:.:.:0,33,495:.
NC_013991.2 3329    .   G   .   18,49   LowQual BaseQRankSum=-1;DP=1818;ExcessHet=3;MLEAC=.;MLEAF=.;MQ=57.63;MQRankSum=0;ReadPosRankSum=1   GT:DP:RGQ   0/0:2:6 0/0:4:12    0/0:2:6 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:1:3 0/0:7:18    0/0:35:99   0/0:13:33

```

It seems like only the "CATT" part of the alternative allele is really inserted...

The mapping was done either with bwa aln or bwa mem, and we've observed this behaviour both when using GATK4 and bcftools for genotyping.

# Answer
> 2 votes
Microsatellites, minisatellites and other type of repetitive variation are a known feature of genomes.

I agree with previous comments in the sense that what you are observing here might as well be real and not just an artefact caused by real repetitive variation in your genome of study, plus read mapping together with variant calling issues at the site.

The latter might be the case, but before getting rid of these alleles I suggest that you inspect the sites and regions around these alleles more carefully: GATK HaplotypeCaller has an option (-- output-bam or similar) that writes a BAM file with the reads as they were used for variant calling after the tool itself performs local re-alignment of reads around the variants. I would use this tool together with the original BAM file in order to inspect the region visually and try to understand it better. For this you can use samtools tview or IGV. The length of your reads becomes somewhat relevant here, but even if you are sequencing around 100 bp, given the small size of these observed alleles you should be able to draw some more solid insights about the variants/errors and how to handle them.

# Answer
> 0 votes
INDELs tend to occur tandemly. Your observation makes them more likely to be real.

---
Tags: vcf, snp, indel
---


# Question
Title: How can I get the coordinates of Hydrogen atoms attached to alpha-carbon?
I am trying to extract the hydrogen atoms attached to the Ca carbon in N no. of proteins. I would be needing the coordinates to calculate its centroid. How should I do it in python ?

# Answer
> 2 votes
OK, tried to figure out a way using Biopython library, I believe there are tens of faster ways, having `1lkk_pdb_test.pdb` as input :

```
ATOM      1  N   LEU A 122      49.672   7.391  11.364  1.00 20.16           N
ANISOU    1  N   LEU A 122     2561   2767   2050    670   -700     42       N
ATOM      2  CA  LEU A 122      50.030   7.561  12.734  1.00 20.25           C
ANISOU    2  CA  LEU A 122     2020   3335   2058   -226   -571    629       C
ATOM      3  C   LEU A 122      48.886   7.480  13.726  1.00 16.62           C
ANISOU    3  C   LEU A 122     1863   2277   1944   -333   -710    706       C
ATOM      4  O   LEU A 122      47.816   8.070  13.517  1.00 15.23           O
ANISOU    4  O   LEU A 122     2132   1556   1887   -122   -845    184       O
ATOM      5  CB  LEU A 122      51.116   8.488  13.074  1.00 21.98           C
ANISOU    5  CB  LEU A 122     1769   3688   2587   -148   -358    491       C
ATOM      6  CG  LEU A 122      51.381   8.978  14.477  1.00 22.50           C
ANISOU    6  CG  LEU A 122     2191   3459   2586   -263   -523    565       C
ATOM      7  CD1 LEU A 122      52.229   8.017  15.285  1.00 23.90           C
ANISOU    7  CD1 LEU A 122     2795   3452   2500   -216   -104    138       C
ATOM      8  CD2 LEU A 122      51.959  10.377  14.502  1.00 26.42           C
ANISOU    8  CD2 LEU A 122     3125   3322   3222    375   -323    398       C
ATOM      9  H1  LEU A 122      48.995   6.786  11.300  1.00 30.24           H
ANISOU    9  H1  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     10  H2  LEU A 122      49.399   8.188  11.020  1.00 30.24           H
ANISOU   10  H2  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     11  H3  LEU A 122      50.392   7.090  10.897  1.00 30.24           H
ANISOU   11  H3  LEU A 122     3829   3829   3829      0      0      0       H
ATOM     12  HA  LEU A 122      50.481   6.686  12.910  1.00 24.30           H
ANISOU   12  HA  LEU A 122     3077   3077   3077      0      0      0       H
ATOM     13  HB2 LEU A 122      51.955   8.066  12.760  1.00 26.37           H
ANISOU   13  HB2 LEU A 122     3339   3339   3339      0      0      0       H
ATOM     14  HB3 LEU A 122      50.982   9.295  12.515  1.00 26.37           H
ANISOU   14  HB3 LEU A 122     3339   3339   3339      0      0      0       H
ATOM     15  HG  LEU A 122      50.493   9.028  14.933  1.00 27.00           H
ANISOU   15  HG  LEU A 122     3419   3419   3419      0      0      0       H
ATOM     16 HD11 LEU A 122      51.771   7.153  15.351  1.00 35.85           H
ANISOU   16 HD11 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     17 HD12 LEU A 122      53.095   7.895  14.844  1.00 35.85           H
ANISOU   17 HD12 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     18 HD13 LEU A 122      52.368   8.382  16.184  1.00 35.85           H
ANISOU   18 HD13 LEU A 122     4540   4540   4540      0      0      0       H
ATOM     19 HD21 LEU A 122      52.806  10.390  14.010  1.00 39.63           H
ANISOU   19 HD21 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     20 HD22 LEU A 122      51.327  10.999  14.083  1.00 39.63           H
ANISOU   20 HD22 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     21 HD23 LEU A 122      52.116  10.649  15.431  1.00 39.63           H
ANISOU   21 HD23 LEU A 122     5019   5019   5019      0      0      0       H
ATOM     22  N   GLU A 123      49.124   6.750  14.812  1.00 14.09           N
ANISOU   22  N   GLU A 123     1166   1937   2055   -326   -424    468       N
ATOM     23  CA  GLU A 123      48.215   6.514  15.890  1.00 12.65           C
ANISOU   23  CA  GLU A 123     1143   1584   1903   -273   -331    -29       C
ATOM     24  C   GLU A 123      46.790   6.330  15.399  1.00 10.88           C
ANISOU   24  C   GLU A 123     1138   1238   1605   -155   -376    -53       C
ATOM     25  O   GLU A 123      45.875   7.059  15.809  1.00 10.93           O
ANISOU   25  O   GLU A 123     1300   1048   1652     37   -173   -144       O
ATOM     26  CB  GLU A 123      48.256   7.537  17.028  1.00 14.54           C
ANISOU   26  CB  GLU A 123     1615   1778   1930   -299   -450   -280       C
ATOM     27  CG  GLU A 123      47.696   6.998  18.341  1.00 19.07           C
ANISOU   27  CG  GLU A 123     2450   2721   1808   -148   -132   -221       C
ATOM     28  CD  GLU A 123      47.620   7.948  19.490  1.00 21.48           C
ANISOU   28  CD  GLU A 123     3199   2805   1858    113    106      6       C
ATOM     29  OE1 GLU A 123      47.461   9.187  19.266  1.00 25.80           O
ANISOU   29  OE1 GLU A 123     4100   2736   2608   1002   -257     -4       O
ATOM     30  OE2 GLU A 123      47.491   7.497  20.672  1.00 22.45           O
ANISOU   30  OE2 GLU A 123     3024   3373   1819    545    690    367       O
ATOM     31  H   GLU A 123      49.919   6.375  14.873  1.00 16.91           H
ANISOU   31  H   GLU A 123     2141   2141   2141      0      0      0       H
ATOM     32  HA  GLU A 123      48.484   5.640  16.294  1.00 15.18           H
ANISOU   32  HA  GLU A 123     1922   1922   1922      0      0      0       H
ATOM     33  HB2 GLU A 123      49.194   7.821  17.171  1.00 17.45           H
ANISOU   33  HB2 GLU A 123     2210   2210   2210      0      0      0       H
ATOM     34  HB3 GLU A 123      47.737   8.336  16.761  1.00 17.45           H
ANISOU   34  HB3 GLU A 123     2210   2210   2210      0      0      0       H
ATOM     35  HG2 GLU A 123      46.784   6.654  18.167  1.00 22.88           H
ANISOU   35  HG2 GLU A 123     2897   2897   2897      0      0      0       H
ATOM     36  HG3 GLU A 123      48.252   6.227  18.618  1.00 22.88           H
ANISOU   36  HG3 GLU A 123     2897   2897   2897      0      0      0       H
ATOM     37  N   PRO A 124      46.576   5.340  14.522  1.00  9.86           N
ANISOU   37  N   PRO A 124      893   1140   1575     13   -285      8       N
ATOM     38  CA  PRO A 124      45.182   5.138  14.041  1.00  8.89           C
ANISOU   38  CA  PRO A 124      875   1151   1229     -6   -227     93       C
ATOM     39  C   PRO A 124      44.244   4.845  15.208  1.00  8.67           C
ANISOU   39  C   PRO A 124      858   1071   1245    -81   -189    -53       C
ATOM     40  O   PRO A 124      44.593   4.143  16.165  1.00  8.98           O
ANISOU   40  O   PRO A 124     1088   1065   1134    -47   -410    127       O
ATOM     41  CB  PRO A 124      45.309   3.946  13.121  1.00  9.41           C
ANISOU   41  CB  PRO A 124     1048   1085   1312     99     -1    -37       C
ATOM     42  CG  PRO A 124      46.555   3.221  13.556  1.00  9.55           C
ANISOU   42  CG  PRO A 124     1054   1161   1282     62   -243      6       C
ATOM     43  CD  PRO A 124      47.489   4.326  13.988  1.00 10.06           C
ANISOU   43  CD  PRO A 124      911   1179   1591    122   -319   -125       C
ATOM     44  HA  PRO A 124      44.868   5.939  13.533  1.00 10.67           H
ANISOU   44  HA  PRO A 124     1351   1351   1351      0      0      0       H
ATOM     45  HB2 PRO A 124      44.517   3.357  13.204  1.00 11.29           H
ANISOU   45  HB2 PRO A 124     1429   1429   1429      0      0      0       H
ATOM     46  HB3 PRO A 124      45.388   4.241  12.180  1.00 11.29           H
ANISOU   46  HB3 PRO A 124     1429   1429   1429      0      0      0       H
ATOM     47  HG2 PRO A 124      46.363   2.605  14.307  1.00 11.46           H
ANISOU   47  HG2 PRO A 124     1451   1451   1451      0      0      0       H
ATOM     48  HG3 PRO A 124      46.943   2.702  12.807  1.00 11.46           H
ANISOU   48  HG3 PRO A 124     1451   1451   1451      0      0      0       H
ATOM     49  HD2 PRO A 124      48.119   4.010  14.683  1.00 12.07           H
ANISOU   49  HD2 PRO A 124     1528   1528   1528      0      0      0       H
ATOM     50  HD3 PRO A 124      48.005   4.678  13.220  1.00 12.07           H
ANISOU   50  HD3 PRO A 124     1528   1528   1528      0      0      0       H
END

```

my code is:

```
import Bio

from Bio.PDB import PDBParser

print('imported Biopython module vesion : ' , Bio.__version__)

parser = PDBParser (QUIET = True)

structure = parser.get_structure('test',  '1lkk_pdb_test.pdb')

atoms = structure.get_atoms()

coordz = []

for i in [i for i in atoms]:

    if i.get_id() == 'HA' :

        print(i.get_parent(), i.get_parent().id, i, i.coord)

        coordz.append(i.coord)

print('coordz list : ')

for i in coordz :
    print(i , type(i))

```

output:

```
imported Biopython module vesion :  1.80
imported Biopython module vesion :  1.80
<Residue LEU het=  resseq=122 icode= > (' ', 122, ' ') <Atom HA> [50.481  6.686 12.91 ]
<Residue GLU het=  resseq=123 icode= > (' ', 123, ' ') <Atom HA> [48.484  5.64  16.294]
<Residue PRO het=  resseq=124 icode= > (' ', 124, ' ') <Atom HA> [44.868  5.939 13.533]
coordz list :
[50.481  6.686 12.91 ] <class 'numpy.ndarray'>
[48.484  5.64  16.294] <class 'numpy.ndarray'>
[44.868  5.939 13.533] <class 'numpy.ndarray'>

```

as per the reason `HA` represent Cα hydrogen, please refer to _

> Westbrook, John D et al. “The chemical component dictionary: complete descriptions of constituent molecules in experimentally determined 3D macromolecules in the Protein Data Bank.” Bioinformatics (Oxford, England) vol. 31,8 (2015): 1274-8. doi:10.1093/bioinformatics/btu789

httpswww.ncbi.nlm.nih.govpmcarticlesPMC4393513pdfbtu789.pdf

Chemical Component Dictionary

as explamples pages for Leucine links to :

https://www.ebi.ac.uk/pdbe-srv/pdbechem/chemicalCompound/show/LEU

and http://ligand-expo.rcsb.org/reports/L/LEU/LEU\_D3L3.gif

lot of trouble to find a diagram:

Just loop over all your proteins in a similar way and store the coords as you need/wish, remember coordinates are retrieved as numpy arrays.

# Answer
> 1 votes
I would import the PDB into the molecular dynamics package GROMACS and use the built-in utilities to calculate the centroid (https://manual.gromacs.org/current/onlinehelp/gmx-rms.html). This is much simpler than trying to write your own utility for this. Other pre-written calculations to perform on your PDB files: https://manual.gromacs.org/current/user-guide/cmdline.html

---
Tags: python, proteins, protein-structure, pdb
---


title: Extract single frame with MDAnalysis

text:
# Question
Title: Extract single frame with MDAnalysis
How can I do to extract a single frame as a pdb file from a trajectory traj.xtc file, by using MDAnalysis?

I've got a universe u=mda.Universe('confout.gro','traj.xtc')

Then I selected a residue by doing select\_atoms, and I would like to extract a single frame as a pdb file

# Answer
> 1 votes
Check out gmx\_trjconv - this command is in GROMACS and accepts several options where you can specify the specific single frame that you want to output/write - out to the new pdb/gro/xtc file. https://manual.gromacs.org/current/onlinehelp/gmx-trjconv.html

---
Tags: python, pdb, gromacs
---


# Question
Title: plot the relation between RMSD and free energy
how can i plot the relation between rmsd and free energy like the image i attached as i need to find the most stable conformation structure for my cycle peptide after finishing REMD simulation

# Answer
You would need to calculate the free energy along a trajectory and the rmsd along that trajectory and then cross-correlate/eliminate variables. There may be a utility which does this but I am not ware of it. https://manual.gromacs.org/current/user-guide/cmdline.html There are two utilities of interest for doing just that. See gmx\_energy and gmx-rmsd. Then you can plot one against the other by elimination of the time/space variables. Check out other calculations you can do with GROMACS: https://manual.gromacs.org/current/user-guide/cmdline.html.
> 0 votes

---
Tags: proteins, gromacs, peptide-shaker
---


# Question
Title: What is the typical host-to-bug DNA ratio found in nanopore sequencing without amplification?
I'm interested in sequencing a human sputum sample using an ONT MinION without performing any type of whole genome DNA amplification or targeted PCR.

Has anyone found a good reference (or anecdotal evidence) for a range of human to pathogen DNA ratios typically found in sputum samples?

Does this depend on whether the patient is infected by a virus or bacteria, and the species that is infecting the patient?

I am seeking a ballpark estimate. Can I expect to find 10x as much human DNA as pathogen DNA, or 100x-1000x?

# Answer
My recollection of a conversation I had with Justin O'Grady is that it's about 99% host DNA for unfiltered samples. There are a few ways to reduce that portion prior to sequencing. He has written a paper that mentions a saponin-based host depletion method that is quite effective:

https://doi.org/10.1038/s41587-019-0156-5
> 1 votes

# Answer
Anecdotally, I have been told that it varies from person to person, also might be affected by disease or hygiene in a particular environment. Here is an article dealing with that problem:

https://journals.plos.org/plosone/article/comments?id=10.1371/journal.pone.0076096

Not sure if they have references useful for you in their paper, but they describe a way of enriching microbial or human DNA. It might be helpful.
> 0 votes

---
Tags: nanopore, covid-19, bacteria, virus
---


# Question
Title: Is there a way I can combine geological data with population genetics data?
So I have ran SAMOVA and made calculations regarding haplotype diversity and nucleotide diversity of several populations. I have also made raster data regarding geographical data like slope, ruggedness etc. Is there a program/package that allows me to combine these to types of data into a map, so I can study if there is a correlation between these indexes and the diversity?

# Answer
> 1 votes
**Update** I looked into modern use of Mantel's test and it will not accept categorical variables. I now remember the null hypothesis of Mantel's ... the association is random between the matrices, i.e. no genetic structure. Therefore *p \< 0.05* excluding degrees of freedom (so in practice its smaller).

The options are:

1. *Fst* implementing that is easy but interpreting is complicated because it's about what constitutes genetic structure

2. General linear modelling;
3. Machine learning for example random forest, SVM etc ...

I would use points 1 and point 3, 2 is old school now. I would use point 1 to datamine the ML.

The key to Fst is that 0.2 comprises population structure for complicated reasons. Implementing ML is not trivial but is the trendiest analysis going. It is not simply implementation its checking and verifying the associations and avoiding overtraining, which would be the major criticism.

The ML would permits the genetic structure to be predicted across broad geographic regions and those predictions are tested and used to reformulate the model. However, you'd need significant expertise in ML (and coding).

---

The technical word for this is **allopatry** and if that doesn't mean much, please do read on the subject because that is central to the theory behind the calculations you are performing. There is a HUGE bank of theory on the subject and is the central alternative to Darwinian selection, yep it's that big. So this is not a passing calculation but a fundamental calculation at the heart of the model of speciation of the organisms within your study. What you want to show is your species does not demonstrate allopatry but niche selection. Ultimately this may not be precisely decipherable (depends on the geographic distances).

Very briefly the analysis to perform now is the **Mantel's test**. This should be within the Arlequin package. Mantel's test will assess if there is a 'cline' in your data, i.e. an association with gene frequency and geographic distance - in an allele frequency package (genetic distance of point mutation otherwise).

Geological data versus geographic distance ... Mantel's test will cover both calculations, however I think that a simple numerical code will be needed for categorical data. I'd need to check that carefully, because I get that wrong the calculation with be wrong. Basically you need to tell the matrix correlation test that this categorical data being tested.

The next level of analysis in allele frequency calculations involves plotting the pairwise (not global) *Fst* against against geographic distance and geological type - if the association in Mantel's test is not linear. If the association is linear the Mantel's will be positive, i.e. null hypothesis will be accept, *p* \> 0.05.

I'll check the location of the calculations in various packages later (not necessarily today). I need to check the categorical variable for Mantel's - I've never done that calculation so am not sure.

**What you are looking for is Mantel's to reject geographic distance and accept geological data type. That would be a decisive and very publishable result.** *Fst* will be the fall-back to shift between the hypotheses.

If you have loads of SNPs there are phylogenetic tree solutions based on modelling point mutation. The next layer of biophylo-geographic analysis involves packages within the Beast framework. I'm sure that would work here because the number of mutations would be low, also this is really just geographic distance rather than niche. If you've loads of SNPs this would be a good analysis to perform however, but for intraspecifies mtDNA particularly involving *cox*1 I personally wouldn't think it would be informative. It is perfectly doable though.

---

A positive Mantel's test is a rare result, but generally allopatry is very common so it is a VERY strong null hypothesis for the data. What I am trying to say is you can't just reject allopatry without testing it. ... well you can but not if you want to publish in e.g. Molecular Ecology.

Good luck with the study and again I'll update the post later. It looks really good and is more than interesting.

---
Tags: phylogenetics, geo, populations
---


# Question
Title: Searching motifs in sequence and their frequencies
This is a two part question. I am searching for a motif, and in that search I wanted to also find the total number of sequences in my FASTA file, but the code I wrote is not yielding that please see the attempt below. The second part is when I have found the matching motifs, as well as the frequencies of the matches how do i go about to exclude the reference sequence matches and only get matches and frequncies from sequence2-sequence5.

```
> Ref
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence1
MAEVRKFTKRLSKPGTAAELRQSVSEAVRGSVVLEKAKLVEPLDYENVITQRKTQIYSDP
LRDLLMFPMEDISISVIGRQRRTVQSTVPEDAEKRAQSLFVKECIKTYSTDWHVVNYKYE
DFSGDFRMLPCKSLRPEKIPNHVFEIDEDCEKDEDSSSLCSQKGGVIKQGWLHKANVNST
ITVTMKVFKRRYFYLTQLPDGSYILNSYKDEKNSKESKGCIYLDACIDVVQCPKMRRHAF
ELKMLDKYSHYLAAETEQEMEEWLIMLKKIIQINTDSLVQEKKDTVEAIQEEETSSQGKA
ENIMASLERSMHPELMKYGRETEQLNKLSRGDGRQNLFSFDSEVQRLDFSGIEPDVKPFE
EKCNKRFMVNCHDLTFNILGHIGDNAKGPPTNVEPFFINLALFDVKNNCKISADFHVDLN
PPSVREMLWGTSTQLSNDGNAKGFSPESLIHGIAESQLCYIKQGIFSVTNPHPEIFLVVR

>Sequence2
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

> Sequence3
MPPRRSIVEVKVLDVQKRRVPNKHYVYIIRVTWSSGATEAIYRRYSKFFDLQMQMLDKFP
MEGGQKDPKQRIIPFLPGKILFRRSHIRDVAVKRLIPIDEYCKALIQLPPYISQCDEVLQ
FFETRPEDLNPPKEEHIGKKKSGNDPTSVDPMVLEQYVVVADYQKQESSEISLSVGQVVD
IIEKNESGWWFVSTAEEQGWVPATCLEGQDGVQDEFSLQPEEEEKYTVIYPYTARDQDEM
NLERGAVVEVVQKNLEGWWKIRYQGKEGWAPASYLKKNSGEPLPPKLGPSSPAHSGALDL
DGVSRHQNAMGREKELLNNQRDGRFEGRLVPDGDVKQRSPKMRQRPPPRRDMTIPRGLNL

>Sequence5
GDDSEWLKLPVDQKCEHKLWKARLSGYEEALKIFQKIKDEKSPEWSKYLGLIKKFVTDS
NAVVQLKGLEAALVYVENAHVAGKTTGEVVSGVVSKVFNQPKAKAKELGIEICLMYVEIE
KGESVQEELLKGLDNKNPKIIVACIETLRKALSEFGSKIISLKPIIKVLPKLFESRDKAV
RDEAKLFAIEIYRWNRDAVKHTLQNINSVQLKELEEEWVKLPTGAPKPSRFLRSQQELEA
KLEQQQSAGGDAEGGGDDGDEVPQVDAYELLDAVEILSKLPKDFYDKIEAKKWQERKEAL
EAVEVLVKNPKLEAGDYADLVKALKKVVGKDTNVMLVALAAKCLTGLAVGLRKKFGQYAG
HVVPTILEKFKEKKPQVVQALQEAIDAIFLTTTLQNISEDVLAVMDNKNPTIKQQTSLFI
ARSFRHCTSSTLPKSLLKPFCAALLKHINDSAPEVRDAAFEALGTALKVVGEKSVNPFLA

```

The code to search the peptides is below

```
import re
import collections
m_list =[]
from Bio import SeqIO
input_file = 'Sequences2.fasta'
num_sequences = 0
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
    num_sequence=+1
    name, sequence = fasta.id, str(fasta.seq)
    matches=re.finditer(r"P[A-Z]{2}P..",sequence)
    for m in matches:
        match_id = f'{m.start()}_{m.end()}_{m.group()}'
        m_list.append(match_id)
        print(name, m.start(), m.end(), m.group())
print("total sequences",num_sequences)
print("frequency",collections.Counter(m_list))

```

It yields the following output below, which I am happy with but now it doesn't give me the total number of sequences:

```
Ref 73 79 PFLPGK
Ref 281 287 PLPPKL
Ref 288 294 PSSPAH
Sequence3 73 79 PFLPGK
Sequence3 281 287 PLPPKL
Sequence3 288 294 PSSPAH
total sequences 0
frequency Counter({'73_79_PFLPGK': 2, '281_287_PLPPKL': 2, '288_294_PSSPAH': 2})

```

The **total sequences** should be **5**. I used the code below to find out. I do not know where I went wrong since it is not giving me that on the more compact code above.

```
sequence = open('Sequences2.fasta')
x = 0
for line in sequence:
    if line.startswith('>'):
        x += 1
sequence.close()
print(x)

```

The second part of my question is how do I modify my code such that, the matches and frequencies of the reference sequence(Ref) are excluded. The idea would be to get an output such as the one below. I do not even know whether that is possible (A disclaimer this second part is a question on a tutorial that I don't know how to do. I am still learning)

```
Sequence3 73 79 PFLPGK
Sequence3 281 287 PLPPKL
Sequence3 288 294 PSSPAH
total sequences 5  #this would be the total number of sequences in file including the refrence
frequency Counter({'73_79_PFLPGK': 1, '281_287_PLPPKL': 1, '288_294_PSSPAH': 1}) #frequency of matches, excluding the reference

```

# Answer
> 3 votes
The first part appears to be a simple bug if there is no copy/paste error in the code:

Original code,

```
num_sequence=0
fasta_sequences = SeqIO.parse(open(input_file),'fasta')
for fasta in fasta_sequences:
    num_sequence=+1 # < its here
print (num_sequence)

```

The code I'd recommend is

```
num_sequence=0
fasta_sequences = SeqIO.parse(input_file,'fasta')
for fasta in fasta_sequences:
    num_sequence += 1 # corrected
print (num_sequence)

```

I will look at the second part later.

# Answer
> 4 votes
for the second part re-arrange the `for m in matches` block to :

```
for m in matches:

        match_id = f'{m.start()}_{m.end()}_{m.group()}'

        print(name, m.start(), m.end(), m.group())

        if num_sequences > 1:       #### added if clause so that appends only after first sequence to m_list

                m_list.append(match_id)

```

---
Tags: python, proteins, sequence-analysis, motifs, multi-fasta
---