# Question
Title: Why are Minimap2 alignments different with CIGAR generation flag?
I am using Minimap2 (v2.26-r1175) in Linux to generate a sequence alignment between the *Streptomyces coelicolor* A3(2) chromosome (`ref.fa`) and the *Mycobacterium tuberculosis* chromosome (`query.fa`). My desired output is a PAF (Pairwise mApping Format) file.

The general way to align reference and query sequences with Minimap2 is the following:

```
minimap2 ref.fa query.fa > approx_mapping.paf
# OR equivalently...
minimap2 -k15 -w10 ref.fa query.fa > approx_mapping.paf

```

You can get Minimap2 to generate custom CIGAR tags (`cg:`) in the PAF by adding the `-c` flag:

```
minimap2 -c ref.fa query.fa > approx_mapping_cigar.paf
# OR equivalently...
minimap2 -c -k15 -w10 ref.fa query.fa > approx_mapping_cigar.paf

```

---

Moving to an R environment \[with the intention of downstream data visualization\], you can quickly assess the number of observations in each object (a.k.a. rows in each dataframe) with `dim()`:

```
file1 <- "approx_mapping.paf"
paf_basic <- read.table(file1, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file1, sep = "\t")))))
dim(paf_basic)
[1] 205  18

file2 <- "approx_mapping_cigar.paf"
paf_cigar <- read.table(file2, sep = "\t", fill = TRUE,
                        col.names = paste0("V",seq_len(max(count.fields(file2, sep = "\t")))))
dim(paf_cigar)
[1] 200  24

```

*Note the unqiue way we have to read the PAF files into R using `read.table()` to ensure no columns are cut off and added to new rows. This is because `read.table()` is designed to examine only the first five lines of an input file to determine the number of columns in the dataframe. Any columns on lines beyond the fifth line that exceed the initially determined count are appended to a new row. In the case of PAF files generated by Minimap2 **with CIGAR tags**, the number of columns can vary per line, which makes this necessary.*

With the **default settings**, I get **205 observations** and with the **addition of the CIGAR flag**, I get **200 observations**.

**Why do I get fewer computed alignments when the CIGAR tags are added?** I would expect the alignments between the reference and query to be identical between these runs.

---

Comparing the PAF files (directly generated by `minimap2`) after sorting (`sort -k3,3 -k4,4`) using `vimdiff` reveals that the alignments computed between the default and CIGAR setting runs are largely *similar*, but not identical:

*Note that the image above only shows the first 44 rows of both files.*

Due to the small differences seen in alignment coordinates between these two dataframes (see fields `$3` and `$4` in the `vimdiff` screenshot), it is **NOT** trivial to determine which five observations (rows) are new alignments computed when using the CIGAR flag `-c`.

---

Plotting these dataframes confirms that the `paf_basic` and `paf_cigar` objects do, in fact, contain some different sequence alignments from each other:

```
library("ggplot2")
ggplot(data = paf) +
  geom_segment(aes(x = V8, xend = V9, y = V3, yend = V4), size = 2) +
  labs(title = "paf_name", x = "Reference (Mb)", y = "Query (Mb)") +
  coord_fixed(ratio = 1) +
  theme_bw()

```

So, it appears as though `minimap2 -c` produces a slightly different alignment than `minimap2`, which is surprising to me and begs the question: **Which alignment is better or preferable?**

# Answer
> 5 votes
The program's author pointed me to the Minimap2 FAQ:

> **1. Alignment different with option `-a` or `-c?`**
>
> Without `-a`, `-c` or `--cs`, minimap2 only finds *approximate* mapping locations without detailed base alignment. In particular, the start and end positions of the alignment are impricise. With one of those options, minimap2 will perform base alignment, which is generally more accurate but is much slower.

---

So, to specifically address the main questions asked here:

1. ***Why are there fewer and non-identical computed alignments when the CIGAR tags are added?***

`minimap2` alignment accuracy differs when running `minimap2 -k15 -w10` *vs.* when running `minimap2 -c -k15 -w10`. It is actually *not unexpected* that the alignment coordinates in each PAF file are only *similar*, not identical, and an increase or decrease in *the number of alignments* between the outputs is simply an artefact of this accuracy discrepancy.

2. ***Which alignment (default vs. with CIGAR tags) is better or preferable?***

Although both alignments report very similar sequence relationships, adding the `-c` flag for CIGAR tag generation is preferable to get more accurate base alignments.

# Answer
> 0 votes
Note also that the "score" can vary dramatically depending on whether -c is used or not. I made a short post about this here https://gist.github.com/cmdcolin/dbb1c6ae291d0e8282cfd0455f908fea it is also alluded to briefly in the PAF.md where score is a required field but "If alignment is not available, column 10 and 11 are still required but may be highly inaccurate."

---
Tags: sequence-alignment, minimap2
---


# Question
Title: Remove sequences from a fasta file with IDs from a text file using Python
a python beginner here. I have a fasta file with 2500+ sequences, and after doing some analysis I want to remove around 200+ sequences based on the matching IDs. Now, I have one fasta file (as sample.fa) and a text file with a list of IDs for the sequences that I want to remove from this fasta file. Below is an example of the fasta sequences. For simplicity, I have shortened the sequences and their respective IDs. The ID in the text file (for the sequences to remove, and in the fast file will be exactly same). The sequences were originally downloaded form NCBI and had very long headers so I have them shortened as at this point I don't need all that info, usually just the sequence ID in the headers is enough.

Sample fasta sequences in fasta file:

```
>WP_017417394.1-Bacillus_amyloliquefaciens_gp.
MRGKKVWISLLFALALIFTMAFGSTSPAQAAGKSNGEKKYIVGFKQTMSTMSAAKKKDVISEKGGKVQKQFKYVDAASAT
>WP_268291751.1-Bacillus_inaquosorum
MRRKKLWISLLFALTLIFTMAFSNMSAQAAGKSSTEKKYIVGFKETMSAMSSAKKKDVISEKGGKVQKQFKYVNAATATL
>WP_003327717.1-Bacillus_atrophaeus
MRSRKLWIGLLFAFTLVFTMAFSSMSPAQAAGKSSEEKKYIVGFKQTMGAMSAAAKKDVIAENGGKLQKQFKYVDAATAT
>WP_206702175.1-Bacillus_vallismortis
MAFSNMSAQAAGKSSTEKKYIVGFKQTMSVMSSAKKKDVISEKGGKVQKQFKYVNAAAATLDAKAVKELKQDPSVAYVEE
>KAF1680298.1-Bacillus_sp.
MKSKKLWISLLFALTLVFTMAFSNMSAQASGKNSTEKKYIVGFKQTMSTAKKKDIISEKGGKVQKQFKYVNAASATLDQK
>WP_172293960.1-Bacillus_sp.
MGKKRLWYSLLTVLLLVFSMAFSHPANAVKSSAKDVKKDYIVGFKSSLKTTSAKKDVIKENGGKVDKHFRSINAAKVTLD
>WP_082194748.1-Bacillus_xiamenensis
MCVKKKNVMTSVLLAVPLLFSAGFGGSMANAETASKSESEKSYIVGFKASATTNSSKKQAVTQNGGTLEKQYRLINAAQV
>WP_048353666.1-Bacillus_glycinifermentans
MKKRSLWLSVLTALLLVFTMAFSNPASAEQPAKDVEKDYIVGFKSSVKTASVKKDVIKESGGKVDKQFKIINAAKATLDQ
>WP_206702858.1-Bacillus_sonorensis
MTALLLVLSTAFSSPASAAQPAKDVEKDYIVGFKSSVKTAAVKKDVIKENGGKVDKQFKIINAAKATLDQEEVKALKKD

```

The list of IDs in a separate text file that need to be removed from the above fasta file:

```
>WP_017417394.1-Bacillus_amyloliquefaciens_gp.
>WP_003327717.1-Bacillus_atrophaeus
>KAF1680298.1-Bacillus_sp.
>WP_048353666.1-Bacillus_glycinifermentans
>WP_206702858.1-Bacillus_sonorensis

```

I searched and tried a few scripts shared elsewhere (e.g. from here), but I couldn't make it work with my files. So all the help will be really appreciated. Thanks in advance.

PS. I am using using python usually with its own IDLE on windows 10 PC.

I did one thing with this code. I removed the '\>' character from the headers file and then this code works but it keeps those headers and their respective sequences. I think it works as intended in this code, but I want to remove those from the sequence file not to keep.

Update: as asked below, I copied this code from another post.

```
from Bio import SeqIO
with open('headers.txt', 'r') as fin:
    to_keep = set([])
    for line in fin:
        to_keep.add(line.strip().split(' ')[0])
with open('output.fa', 'w') as fout:
    for record in SeqIO.parse('sample.fa','fasta'):
        if record.id in to_keep:
            SeqIO.write(record, fout, 'fasta')

```

# Answer
After modifying the code above in the original post, this one works the best for me. I got it working after putting 'not' in the second last line instead of "if record.id in to keep" to "if record.id not in to keep":

```
from Bio import SeqIO
with open('headers.txt', 'r') as fin:
    to_remove = set([])
    for line in fin:
        to_remove.add(line.strip().split(' ')[0])
with open('output.fa', 'w') as fout:
    for record in SeqIO.parse('sample.fa','fasta'):
        if record.id in to_remove:
            pass
        else:
            SeqIO.write(record, fout, 'fasta')

```
> 2 votes

# Answer
You can use something like seqkit: https://bioinf.shenwei.me/seqkit/

seqkit has a `grep` sub-command that is one of the most flexible I've ever seen. Your use case is quite simple, you can use `seqkit grep -v` like so:

```
seqkit grep -v -n -f list_of_IDs_to_exclude.txt in.fasta > out.fasta

```
> 5 votes

# Answer
This should work

```
import argparse

ap = argparse.ArgumentParser(description=__doc__, formatter_class=argparse.ArgumentDefaultsHelpFormatter)
ap.add_argument('--fasta', required=True, type=str, help='Fasta file to extract from')
ap.add_argument('--headers', required=True, type=str, help='Headers to remove, one per line.')

conf = ap.parse_args()
f = conf.fasta
h = conf.headers
header_set = set()

with open(h) as header_file:
    for line in header_file:
        line = line.rstrip()
        header_set.add(line)

with open(f) as fasta_file:
    remove_sequence = False
    for line in fasta_file:
        if '>' in line:
            header_id = line.split()[0].replace(">", "")
            if header_id in header_set:
                remove_sequence = True
            else:
                remove_sequence = False
                print(line.rstrip())
        elif not remove_sequence:
            print(line.rstrip())

exit()

```

Save it as python script and then run

`python remove_fasta.py --fasta sample.fa --headers listofheaders.txt > cleaned.fasta`
> 2 votes

# Answer
I would argue writing a whole python script for this is overkill. Assuming you are using a\*nix based operating system like Linux or macOS, you can use the two scripts I have posted here, `FastaToTbl` and `TblToFasta`, and just `grep` out the ones you do not want:

```
FastaToTbl sample.fa | grep -v -xFf listofheaders.txt | TblToFasta > cleaned.fa

```
> 2 votes

---
Tags: python, fasta, sequence
---


# Question
Title: Longstitch error make: command: Command not found *** No rule to make target
I installed Longstitch and ran the test script with no issues. The output files matched the expected output files. But when I am now trying to run Longstitch on my own data I am getting this error.

```
(medaka) operations@LabPC1:~/medaka/longstitch-1.0.5$ ./longstitch run draft=~/medaka/~medaka/medaka_results/consensus.fa reads=~/medaka/combined.fastq G=4.5e8
make: command: Command not found
make: command: Command not found
make: *** No rule to make target'/home/operations/medaka/~medaka/medaka_results/consensus.fa.fa', needed by'/home/operations/medaka/~medaka/medaka_results/consensus.fa.cut250.tigmint.fa'.  Stop.

```

I changed the file name from .fasta to .fa for my draft assembly because it requires a .fa file, but now I am getting ".fa.fa" What is wrong here?

# Answer
> 1 votes
Could you download it from `conda` and try again? Its here, https://docs.conda.io/en/main/miniconda.html

Then,

```
conda create --name myenv
conda activate myenv
conda install -c bioconda longstitch
longstitch ...

```

Conda avoids all the build stuff.

---

Notes

If you find that "conda" is slow and doesn't ever load.

There is an alternative to conda called **mamba**, this is fast and designed to address this exact issue. I have never needed to use it so the following is a rough guide only.

Download and install via

```
conda create --name myenv
conda activate myenv
conda install -c conda-forge mamba

```

The documentation is here, https://mamba.readthedocs.io/en/latest/user\_guide/mamba.html

I've never used it but the idea is ...

```
mamba create -n mambaenv <list of packages>
mamba activate mambaenv

```

The list of packages could be left blank for now, then,

```
mamba install "Python>=3.10.0" numpy pandas pathlib
mamba install longstitch

```

Can also do

```
mamba install -n mambaenv longstitch

```

**NOTE** longstitch exists strictly in -c bioconda. You do not seem to search and install via channels in mamba

You search for a package and its dependencies via

```
mamba repoquery search longstitch
mamba repoquery depends -a longstitch

```

It is important to place a virtual environment between the system Python and package Python (don't want a system corruption).

---

Finally ... For `fa.fa` ... hmmmm, try `ls *.fa` and check for `.fa.fa` files.

# Answer
> 1 votes
The main issue I can see is the filenames, and should be easily solvable.

Simply put, LongStitch is *very* particular about filenames. It expects the input draft to have the extension .fa, and the input reads to have any from a range of extensions (your .fastq is perfectly fine, though).

The main issue is that LongStitch expects the extensions to be *missing* from the command. So if your draft is draft is called `consensus.fa` and your reads are `combined.fastq`, then the command would include something like `draft=consensus reads=combined`, not `draft=consensus.fa reads=combined.fastq`.

(unsure if it accepts input from different directories, or whether you'll need to move your input to the same directory as well)

More generally, when having issues with a particular tool which is available via github, I would encourage you to have a look at the 'Issues' or 'Help' section. What I normally do is to delete the 'is:issue is:open' terms in the search bar, which allows me to search for my bug among older, potentially solved issues (e.g. https://github.com/bcgsc/LongStitch/issues/54), and also generally to check that the software is still supported.

Good luck!

(\*that said, M\_'s suggestion of using mamba is generally a good one, and has helped me with another LongStitch bug)

---
Tags: sequence-analysis, assembly, wgs, long-reads, scaffold
---


# Question
Title: How to annotate a bacterial genome automatically?
We recently got some shotgun sequencing results from a soil bacterium. And we have obtained some contigs which have some genes of interest.

Is there a way to automatically annotate the whole sequence (one single contig of around 25 kb size) and mark the genes on this stretch of DNA? I have never done this before so any help would be appreciated. The ultimate goal would be to know what genes are on this stretch of DNA and see if they work together in a pathway. But the first problem is:

How to annotate a stretch of DNA with many genes on it that would code for a protein/enzyme? so that I get a map of genes on this stretch of DNA. Thanks.

Edit:

This is a soil isolate and we don't have a full genome sequence. Presumably, it is a Rhodococcus species, which might be similar. From the one gene we know, there are no exact homologs in the database so most probably it is a new strain. The strain works symbiotically with another strain which we are also working on and later I want to predict what proteins maybe involved in the catalysis of a specific compound. It is not a one step process so possibly multiple enzymes/proteins are involved. I tried UGENE and UGENE can identify multiple ORFs in this particular stretch, but there are hundreds of ORFS (on both strands) and I find it tedious to blast on NCBI one by one.

I didn't want too many questions in one post, so at the moment just focus on the identification of all the putative genes on this one contig would be good enough to start.

Thanks again.

# Answer
There are a number of annotation pipelines that you could use to get fairly standardized annotation - Prokka has been around for a while (paper published in 2014) but is still used frequently. There's also Bakta, which may be better at catching short ORFs and providing more detailed functional annotation, but requires a fairly large database.
> 5 votes

---
Tags: genome, sequence-annotation
---


# Question
Title: Volcano Plot: How to show unique gene names in the plot?
I'm wondering how can I show unique gene names in the volcano plot.

My data looks like,

```
head(DMCf1_huMetGNc)
     chr    start    end.x       pvalue       qvalue meth.diff        gene
152 chr1 24909805 24909805 1.915400e-03 0.0168087962 -52.00000 RUNX3;RUNX3
168 chr1 24931440 24931440 8.645625e-03 0.0489503201  28.98551 RUNX3;RUNX3
179 chr1 26862096 26862096 1.308307e-04 0.0023224596 -42.48366         SFN
187 chr1 26914178 26914178 4.508416e-06 0.0001786963 -57.97872       NR0B2
188 chr1 26914203 26914203 5.337425e-05 0.0011910452 -47.30424       NR0B2
194 chr1 35573100 35573100 1.433689e-03 0.0135780139 -42.04545      TFAP2E

```

The original data file is quite large, here under the `gene` column some of the rows have multiple same gene names separated by `;`. How can I show only unique gene names.

No pattern of gene names exists when multiple names are present. In certain cases, the same gene is repeated, while in other cases, distinct genes are present in the same cell. The majority of the time, a single gene exists in a cell.

```
library(ggplot2)
library(ggrepel)
    v1<-ggplot(data=DMCf1_huMetGNc, aes(x=meth.diff, y=-log10(pvalue) )) + geom_point()
    v1+geom_text_repel(aes(label= gene))

```

Here, you may see how the FGFR2 gene is repetitive in the plot that I want to remove and the same name from a single row should appear only once in the plot.

# Answer
> 1 votes
Based on your comments, you wish to obtain a semi-colon separated list of unique values in the gene column. You can create a new column containing these using:

```
xyz <- data.frame(chr="chr1", gene = c("RUNX3;RUNX3", "RUNX3;RUNX3", "SFN", "NR0B2", "NR0B2", "TFAP2E", "A;B;A"))

xyz
   chr        gene
1 chr1 RUNX3;RUNX3
2 chr1 RUNX3;RUNX3
3 chr1         SFN
4 chr1       NR0B2
5 chr1       NR0B2
6 chr1      TFAP2E
7 chr1       A;B;A

xyz$gene_uniq <- sapply(xyz$gene, function(x) paste(unique(unlist(strsplit(x, ";"))), collapse=";"))

xyz
   chr        gene gene_uniq
1 chr1 RUNX3;RUNX3     RUNX3
2 chr1 RUNX3;RUNX3     RUNX3
3 chr1         SFN       SFN
4 chr1       NR0B2     NR0B2
5 chr1       NR0B2     NR0B2
6 chr1      TFAP2E    TFAP2E
7 chr1       A;B;A       A;B

```

Or you can use dplyr. I'm adding a column called `entry_num` below since my example has non-unique combinations of chr and gene that we need to retain.

```
library(dplyr)
library(tidyr)

xyz <- data.frame(chr="chr1", gene = c("RUNX3;RUNX3", "RUNX3;RUNX3", "SFN", "NR0B2", "NR0B2", "TFAP2E", "A;B;A"))
xyz %>% mutate(entry_num = 1:nrow(xyz)) %>% separate_rows(gene, sep = ";") %>% unique() %>% group_by(chr, entry_num) %>% summarize(gene = paste(gene, collapse=";")) %>% select(-entry_num) %>% as.data.frame()

   chr   gene
1 chr1  RUNX3
2 chr1  RUNX3
3 chr1    SFN
4 chr1  NR0B2
5 chr1  NR0B2
6 chr1 TFAP2E
7 chr1    A;B

```

---
Tags: r, genomics, ggplot2, methylation
---


title: Classical alignment with HISAT2

text:
# Question
Title: Classical alignment with HISAT2
*This question was also asked on Biostars*

I have this script

```
#!/bin/bash

# Defining the number of threads to be used
threads=20

# Alignment of sequences to the reference genome using 'HISAT2'
echo -e "\nNow Running Classical Alignment-based with HISAT2...\n"

mkdir -p hisat2

# Download the human reference genome
wget -c ftp://ftp.ensembl.org/pub/release-100/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz -O hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz

# Download the human annotation file
wget -c ftp://ftp.ensembl.org/pub/release-100/gtf/homo_sapiens/Homo_sapiens.GRCh38.100.gtf.gz -O hisat2/Homo_sapiens.GRCh38.100.gtf.gz
gunzip hisat2/Homo_sapiens.GRCh38.100.gtf.gz

# HISAT2 indexing to build the reference genome index
hisat2-build -p $threads hisat2/Homo_sapiens.GRCh38.dna.primary_assembly.fa hisat2/Homo_sapiens.GRCh38v3_hisat2.idx

# HISAT2 Alignment
for sample in $(cat sample_ids.txt)
do
    hisat2 -p $threads -x hisat2/Homo_sapiens.GRCh38v3_hisat2.idx \
        -1 trimmed_data/${sample}_1_val_1.fq -2 trimmed_data/${sample}_2_val_2.fq \
        -S hisat2/${sample}_hisat2.sam

    # Compress the SAM files to binary format, sort them, and index them.
    samtools view -Sb hisat2/${sample}_hisat2.sam | samtools sort > hisat2/${sample}_hisat2_sorted.bam
    samtools index hisat2/${sample}_hisat2_sorted.bam
rm hisat2/${sample}_hisat2.sam  # Remove the .sam files for storage purposes
done

# Quantifying Aligned reads using the Subread package's script 'featureCounts'
for sample in $(cat sample_ids.txt)
do
featureCounts -T $threads -a hisat2/Homo_sapiens.GRCh38.100.gtf \
        -o hisat2/${sample}_hisat2_counts.txt \
    hisat2/${sample}_hisat2_sorted.bam
done

# Generating analysis report with MultiQC
echo "Running MultiQC on HISAT2 alignment results ........."
multiqc -o multiQCreports hisat2
echo "MultiQC on HISAT2 alignment results completed successfully."

# Display message to show that the script is complete
echo "Script completed!"

```

However, I have receieved an error below

```
Load annotation file Homo_sapiens.GRCh38.100.gtf ...                       ||
||    Features : 1378275                                                      ||
||    Meta-features : 60683                                                   ||
||    Chromosomes/contigs : 47                                                ||
||                                                                            ||
|| Process BAM file SRR19170514_hisat2_sorted.bam...                          ||
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam
ERROR: Paired-end reads were detected in single-end read library : hisat2/SRR19170514_hisat2_sorted.bam

```

Question: What can I do to solve this?

# Answer
I believe the error is from `featureCounts` and the solution would be to add the `-p` flag to the `featureCounts` call, stating that you are providing paired-end data.

EDIT: You might want to check this `r/bioinformatics` (reddit) thread out.
> 6 votes

---
Tags: rna-seq
---


# Question
Title: Can I interpret 15x coverage area as a part of one strand is missing?
This is a BAM file from a 30x coverage Nebula Genomics test opened with IGV.

I found an area (1:30169255-30169600) where the coverage is less than 15x. Can I interpret it as a part of one strand of Chromosome 1 is missing?

# Answer
## UPDATE:

Seeing a different question where similar vocabulary is used: it would not be one *strand* that would be lost. single stranded DNA is not very happy in the nucleus, and would either lead to something seriously bad or be gap-filled to make it dsDNA again.

It would instead be that one of two chromosomes has a deletion at that position, which would be a classic CNV.

## ORIGINAL ANSWER:

No.

There are a variety of copy number variant (CNV) calling tools that put a lot of effort into statistically robust ways of asking this question, all of which are a lot more complex than just looking at raw coverage. I'd recommend looking more into them.

In terms of explanations for lower coverage in some regions that aren't CNVs: sequenceability of DNA is not even, replication of DNA is not even, mappability of reads is not even, etc. You need to reject those hypotheses before you go to variants.
> 1 votes

---
Tags: wgs, indel, igv, chromosomes
---


# Question
Title: Single Cell RNASeq Knee Plot - Collapsed vs Non-collapsed Reads
I am trying to adapt combinatorial index reads with some modifications for long read sequencing; it's a custom approach, not a commercial kit.

I have a question about using Knee Plot for identifying Barcode white list for non-10X data. While generating the Knee Plot barcode which barcode frequencies should be used: without collapsing PCR replicates to after collapsing PCR replicates?

I have used UMI-Tools (https://umi-tools.readthedocs.io/en/latest/QUICK\_START.html) and Pacbio groupdedup based approach. For the knee plot generation I used both collapsed and non-collapsed reads. Since its long reads and coverage is not great; bumps are not the perfect one but both inputs generate it. I am not sure which input I should consider to remove the backgrounds.

Thanks

**EDIT:** *Can the down-voters please leave a note why they think this not a good query to ask?*

# Answer
> 0 votes
If you have UMIs, then cell/UMI/mapped transcript triplets should be collapsed into unique triplets.

---
Tags: scrnaseq, barcode, umi
---


title: How to visualize FIMO outputs?

text:
# Question
Title: How to visualize FIMO outputs?
I used FIMO to scan several of my promoter sequences for TF binding motifs. The output .tsv file gives me the motif name/sequence plus the start and stop location, how can I visualize these outputs to get a graph that shows each motif as differing length rectangles along each promoter sequence, along with a gradient that shows the p-value?

# Answer
Potential workaround: use MAST.

As far as I understand, 'under the hood', so to speak, FIMO and MAST work pretty much the same. The difference might be just that FIMO focuses on scoring motifs, while MAST is scoring sequences (https://groups.google.com/g/meme-suite/c/qMSj6AdFxoM)

Example of FIMO output: https://web.mit.edu/meme\_v4.11.4/share/doc/examples/fimo\_example\_output\_files/fimo.html

Example of MAST output: https://web.mit.edu/meme\_v4.11.4/share/doc/examples/mast\_example\_output\_files/mast.html

There might be small differences in the output of the tools, but then it's up to you to see whether this is relevant for your data.
> 0 votes

# Answer
If the sequences you are scanning with FIMO are supported by the UCSC or Ensembl genome browsers you can upload the GFF file generated by FIMO as a custom track: UCSC Custom Track, ENSEMBL Custom Track. Note that for the Ensembl website the FIMO GFF file is in GFF3 format.
> 0 votes

---
Tags: bed, motifs, tsv
---


# Question
Title: How can I programmatically add Hydrogen to a PDB structure using BioPython?
I need to programmatically add Hydrogen atoms to a PDB file (say, 16PK.pdb) using BioPython.

What algorithm should I use for this?

# Answer
> 4 votes
Not actually and answer as you cannot do it with `biopython.PDB`. But there are may other tools.

Several are online, for example: https://zhanggroup.org/HAAD/ by the folk who wrote I-Tasser. Nearly every MD pipeline has their own protonation algorithm. The only hydrogens that are an issue are histidine protonation —`HIE` vs. `HID`— and protonation states of catalytic residues (`CYZ` vs. `CYS` etc).

The PyMOL python module is really good and can do a lot. Only issue is that the easiest way to install is via conda.

```
conda install -c conda-forge -c schrodinger pymol-bundle

```

with which one can:

```
import pymol2
pymol2.PyMOL() as pymol:
    pymol.cmd.load('.pdb', 'myprotein') # file
    # pymol.cmd.fetch('', 'myprotein') # PDB code
    pymol.cmd.h_add()
    pymol.cmd.save('.pdb')

```

One can find neighbours to a given atom too (see this Q&A), select by secondary structure, define secondary structure (`pymol.cmd.dss()`) etc. A caveat is that, bar for some commands like `pymol.cmd.get_model(selection).atom`, one does not interact directly with the atoms and everything is done via selections.

However, as this Q follows on from this unanswerable one and in turn from several others about SS. I am guessing it is to determine the hydrogen bond distances between different sheets. These have a backbone amide nitrogen (name `N`) donating a hydrogen (`H`, not `HA` etc.) to an accepting backbone amide oxygen of another sheet (`O`). Unless you want an exact distance, you could simply use the distance between these heavy atoms, which would be slightly larger than the `H-N` bond plus the `H•O` hydrogen bond and the oxygen may not be coplanar with the planar amide.

# Answer
> 0 votes
`conda install -c conda-forge openbabel`

```
from openbabel import pybel as pb
babel_mol = pb.readfile(format="pdb", filename=path_to_pdb).__next__()
babel_mol.addh()
babel_mol.write(format='pdb', filename=path_to_babel_pdbh, overwrite=True)

```

---
Tags: python, proteins
---


# Question
Title: Do additional peaks in percent GC of PacBio gDNA reads indicate contamination?
I have two sets of PacBio reads from genomic DNA of an Aspergillus species that were made from separate preps of the culture. One of them has two additional peaks at 38% and 60% in the percent GC histogram produced by FastQC. Do these additional peaks indicate that there is contamination in that sample?

# Answer
The logic I would us is that if they come from the same genome, coverage peaks will go in multiplications of 2 - 1n, 2n, 4n, ..., while if it's a contamination, the coverage will go wild. So, I would recommend to look at the coverage vs GC to decide where it comes from.

You could either assemble the genome and use something like Blobtools to check coverage, GC and putative origin of sequences (and you should do that anyway).

Perhaps a bit simpler would be to look at it using k-mers, but that assumes you have reads with low error rate (In your case PacBio CCS or HiFi). If you do, you can use KAT that can do GC vs coverage plot - they intended it to look at GC bias in sequencing, but the vary same plot should show you if the three GC peaks have the same coverage or not.
> 1 votes

# Answer
I would suggest plugging your reads into a fast k-mer matching tool like mash screen to get a basic idea of where your reads come from, if you suspect contamination.

An issue with just looking at GC content is that some organelles/plasmids/etc. such as mitochondria have quite divergent GC content, and the different populations you observe might just indicate random fluctuations of copy number of these DNAs due to different sample preps. Not sure that would explain your plot but it does give me pause.

As suggested above, assembly is obviously a good idea if you are ok with the compute cost, and then you can use something like SPAdes or blobtools to pull out contigs from different origins.

KAT and similar tools as suggested are handy for hypothesis free evaluation of different populations in your reads, but they don't directly address the issue of contamination.

And obviously measure the quality of the reads to establish that they are ok and nothing weird happened there, they aren't all Gs or something.
> 1 votes

# Answer
Yes, it looks odd, but you might need to assemble the sequences to tell what is happening. The fact that the peaks are symmetrical around the presumably *Aspergillus* peak around 50% GC is strange. Usually contamination appears as one peak or shoulder on either side of the host DNA peak rather than symmetrical peaks.
> 0 votes

---
Tags: genome, pacbio
---


title: Polishing PacBio or ONT with Illumina

text:
# Question
Title: Polishing PacBio or ONT with Illumina
Which tool would be good to polish PacBio or ONT with Illumina?

Thank you in advance

# Answer
> 1 votes
PacBio and ONT are very different technologies; it would be helpful to be more specific. Regardless, I can only give my understanding as it relates to ONT reads....

TLDR: use medaka. Don't use Illumina reads.

Nanopore reads called using the most recent super-accuracy basecallers have a very high consensus accuracy (q40-q50) for SNPs and non-homopolymeric INDELs after nanopore-only correction. ONT provide their own tool, medaka, which uses assembly-mapped nanopore reads to correct the assembly. Medaka uses the base calling model from Guppy to predict the most likely sequence, given the input reads and the assembly.

If you're using Illumina reads to polish a nanopore assembly, it should *only* be used to correct INDEL errors (and ideally, only homopolymer errors). The two most commonly-used tools that can correct nanopore reads using Illumina reads are Racon and Pilon. Based on the recency of the commit histories, I'd guess that Racon would be the most likely to work well. It looks like Racon can do both assembly correction *and* read correction.

# Answer
> 1 votes
Why don’t you use a hybrid assembly which pretty much does it by default? Tools like OPERA-MS assemble the short reads with MEGAHIT and then uses nanopore to build bigger contgs. Otherwise, I recommend pilon which probably the most popular tool.

# Answer
> 0 votes
For PacBio HiFi reads (or other techs of comparable quality), it is not recommended to polish with Illumina at all. You run the risk of collapsing true heterozygosity/variation in your underlying genome and yielding a haploid assembly that does not represent either haplotype.

If that underlying variation is an issue because you need a haploid assembly, you are free to use a deduplication technique such as purge\_haplotigs or similar (reviewed here) to obtain a primary set of contigs/scaffolds. Assigning p and q haplotypes is something that you can do with a tool like falcon-phase (if you have Hi-C) or a similar technique (I haven't been following the literature for a couple years, possibly there's something a lot better now!).

For a recent review of recommended strategies you can read the T2T consortium benchmarking paper on the topic.

---
Tags: nanopore, illumina, pacbio
---


# Question
Title: align to the whole hg38 genome, then split bam and collect metrics on each bam issue
I am doing some coverage analysis on deep sequencing human genome. I first align to the whole hg38 genome, then split bam to each chrom, and collect metrics on bam for each chom. I first split using samtools. and then tried bamtools.

```
samtools: samtools view -bo ${pair_id}.\\\$chr.bam $bam \\\$chr
bamtools:bamtools split -in $bam -reference

```

When I collect metrics I am using the reference fasta from each chrom. The picard metrics complains the index dictionary size is different: Sequence dictionaries are not the same size (25, 1). The picard method I used is: MarkDuplicates, CollectAlignmentSummaryMetrics, CollectGcBiasMetrics, CollectInsertSizeMetrics, CollectWgsMetrics.

My data is really huge, 70G for a sample, that is why I would like to split when doing the picard metrics. Seems picard metrics collection use a lot of memory. So for the 70G human fastq(30X), how big the memory is needed for running picard wgs metrics?

How can I deal with this? Can I modify the split bam header? How to do it? Or not split it into chromosmes, rather, using an interval list when collecting the metrics? will that use less memory?

# Answer
> 1 votes
This can be done, but you have to split the BAMs while keeping the full header (with the dictionary) in place for all the shards. This can be done with gatk PrintReads. One thing to remember is that you all need to have one shard with all the "unmapped" reads. Otherwise you will not have the complete information. You can get the unmapped reads by putting `unmapped` into the `-L` argument.

That said, MarkDuplicates cannot be run on shards since it needs both reads from every pair and some of the read pairs might be split. So you'll need to mark-duplicates on the entire 70GB bam (sorry).

You could split the BAM into shards but only the "proper pairs" ( filter after GATK with `samtools view`, but only on the sam flag so you don't mess with the dictionary) which will make sure that your pairs are on the same contig. Then put all the not proper pairs into their own shard.

You can also do the proper/improper filtering with GATK, but it involves read filters which are a bit more complex than samtools flag filters.

---
Tags: sequence-alignment, picard
---


# Question
Title: probability of finding a 5 amino acids in a row within a proteome
How to calculate the probability of finding two proteins that share a 5 amino acid long motif from a proteome of around 1067 proteins that have an average length of 65 residues. The probability of a specific sequence that is completely the same is 1067/20^65

# Answer
> 2 votes
Its just significant at 2%

> 0.055 . 65 . 1067 = 0.021 (*its actually 0.02 - see note*)

Where critical probability is 0.05, i.e. 5%

There are 20 aa so 1/20 probability of random occurrence. The oligopeptide is 5 amino acids long, hence to the power 5. Obviously there are 1067 proteins and the probability is per protein, thus if these were independent its a straight multiplication because that is how independent events are measured.

**Note** The 65 is really 60 because it's a sliding window moving 1 amino acid at a time and the proteins are only 65 in length. The actual value is 60 rather than 65 because at the 3' terminal of the sequence the sliding window can't capture the remaining 4 amino acids, because once the oligopeptide \< 5 the calculation is no longer valid. I used 65 so you could see where the value is originating.

Obviously the proteins in question are not independent so you need to state the assumption, i.e. assuming independence.

---
Tags: sequence-alignment, proteins, statistics, sequence-analysis, homework
---


# Question
Title: what percentage of the human genome is MAPQ=0?
When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, what percentage of the non-N fraction of the human genome is MAPQ=0? This is, what part corresponds to regions that can't be uniquely mapped with 2x150bp reads.

And, how many genes are affected by MAPQ=0 regions?

I presume the numbers will dance around depending on fragment sizes, read quality, etc. but I am happy with some starting numbers.

# Answer
Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper is not necessarily MAPQ=0 to another.

That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a useful answer. For 150bp reads, each reference position is covered by 150 reads. What if 50 of them have no other exact hits elsewhere, but the rest 100 have? Is this position a repeat or not? In addition, what if the 50 are unique only because one mismatch? If there is a variant at that mismatch, the 50 would become repeats or mapped elsewhere.

My preference is to say a position is "unique" under k-long reads if over k/2 reads overlapping the position have no other perfect or 1-mismatch/1-gap hits elsewhere in the genome. Under this definition, 79.3% of human genome are unique for 35bp reads. 92.4% for 75bp reads. I don't have the number for 150bp reads. I guess it will be around 95%. Empirically, 94-95% of human genome is callable with 100bp paired-end reads.

As to other measurements, the most common one is the fraction of reads that has an exact hit elsewhere. You can use Fabio's method. It gives a good enough estimate once you simulate over 1 million reads. This fraction is around 86% for 35bp reads and 95% for 75bp reads, as I remember. The problem with this approach is this fraction is not very informative to variant calling due to the issues I talked about. Another way is to use RepeatMasker. It is worse. RepeatMasker masks 50% of human genome, but excludes segmental duplications where short reads can't be confidently mapped.
> 6 votes

# Answer
As has been said before, mappability to the 'human genome' depends on a number of factors, among these the reference version and type of reads, for which you are interested in GRCh38 and 2x150bp reads. Although I am not aware of numbers accounting for these particular reference and type of Illumina reads, the 1000 genomes project has provided the community with a similar and close estimation that you might be interested in considering regarding your inquiry.

Similar to your question, the 1K project estimated 'the proportion of the human genome that is less accessible to short reads'. In these estimates the human genome is GRCh37 and the types of reads in question are mostly 2x Illumina with a mixture of lengths with the longest being up to 250bp. In these estimates each base in the human genome is considered (and marked) 'less accesible' according to these criteria:

L - depth of coverage is much lower than average

H - depth of coverage is much higher than average

Z - too many reads with zero mapping quality overlap this position

Q - the average mapping quality at the position is too low

Each of these criteria has "standard" and "strict" thresholds for a base to be considered - or not - in each category. You can read more in the link below:

ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/accessible\_genome\_masks/README.accessible\_genome\_mask.20140520

According to the strict thresholds, the human genome has about 16.8% of "Z" and 3.1% of "Q" bases, respectively. Considering the "Z" and "Q" criteria as a proxy for ~ mapq=0, about 19.9% of the human genome can not be uniquely mapped.
> 5 votes

# Answer
I would suggest you to run some simulation. You can use wgsim to simulate reads that are as much similar to what you want an answer for. I don't think there are faster methods.
> 1 votes

# Answer
You might want to check out this paper:

They have collected a list of "dark" (regions that have few reads aligning to them) and "camouflaged" (regions that have ambiguous alignment) regions of the human genome:

> Here, we identify regions with few mappable reads that we call dark by depth, and others that have ambiguous alignment, called camouflaged.
>
> \[. . .\]
>
> Other dark regions arise, not because the sequencing is inherently problematic, but because of bioinformatic challenges. Specifically, many dark regions arise from duplicated genomic regions, where confidently aligning short reads to a unique location is not possible; we term these regions as “camouflaged”. These camouflaged regions are generally either large contiguous tandem repeats (e.g., centromeres, telomeres, and other short tandem repeats), or a larger specific DNA region that has been duplicated (e.g., a gene duplication) either in tandem or in a more distal genome region. In fact, many genes in the human genome were duplicated over evolutionary time and are still transcriptionally and translationally active (e.g., heat-shock proteins) \[3–9\], while others have been duplicated, but are considered inactive (i.e., pseudogenes). Regardless of whether the duplication is active, however, any genomic region that has been nearly identically duplicated and is large enough to prevent sequencing reads from aligning unambiguously will be “dark”, because the aligner cannot determine which genomic region the read originated from.

The authors have made their scripts available on github, along with bed files for the relevant regions. The answer to your question, i.e. a specific percentage, is hard to give since it varies enormously with the genome build and sequencing technology, but this paper should have all the details you want. Have a look at Table 1:
> 1 votes

---
Tags: genome, read-mapping, human, illumina, mapq
---


title: Process gets killed before completing

text:
# Question
Title: Process gets killed before completing
I have a large tsv file (~15 GB) that is a kmer count matrix and I am trying to calculate jaccard distance of the counts but every time I try to run the following code the process ends before it's done with a `killed` error message in the STDOUT:

```
jaccard = pd.DataFrame(squareform(pdist(kmer_matrix, 'jaccard')))
jaccard.to_csv(matrix_path+kmer+"_"+method+"_jacDistance_noFilter.tsv", sep='\t')

```

This process works really well with smaller tsv files so I know the code works but I wanted to know if anyone knew what would help run the process quicker/completely. I know the problem must be with my RAM/CPU but is there a way to maybe break up the matrix to calculate everything in batches to make it easier to work with?

# Answer
You've overloaded the memory. One route is to deploy "chunking". pandas will do this. This is complicated to implement code wise however.

My suggested solution (because it simpler) would be to do this via `pandas2` which uses Apache Arrow via PyArrow. You need to specify which method is being used at the onset of declaring the dataframe. Apache Arrow is for lifting large amounts of data. If this works there's no complicated coding you just `import pandas2 as pd` and in the dataframe construction describe the precise heavy lift method (I think there's an additional approach besides PyArrow).
> 3 votes

---
Tags: k-mer, errors, distance-matrix, troubleshooting
---


title: pacbio demultiplexing issue

text:
# Question
Title: pacbio demultiplexing issue
I was wondering if there is a way to provide pacbio lima demultiplexing only a barcode, and let lima demux sample even if only one end has the barcode, not both?

# Answer
PacBio lima need both 5'end and 3' end barcodes. You can have a unique barcode in one end and multiple on the other end; but it need two barcodes.
> 0 votes

---
Tags: pacbio, demultiplexing
---


# Question
Title: Fastq file with very high per base sequence quality and no box-whisker plot
I have my data from NCBI. I am having difficulty making sense of and correlating my fastq file read 1 and read 2 and the fastQC report. As shown below, I have two paired-end read results each displaying 127 bp in length; however, my results on the fastQC report of per base sequence quality show 91 and it does not show the interquartile range (box-whisker plot). What could be the problem?

Below is an excerpt display of my fastq file and then the fastQC report:

```
@SRR20708369.1 NS500749:581:HGV7VBGXB:1:11101:17225:1048 length=127
AGCTATCACAGCANGGTACTTCCCAGCGTTACGGGANTGGCCGTACGCCTAACCGCTAACATTACTGCAGGCAACCTACTCATGCACCTAATTGGAAGCGCCACCCTAGCAATATCCACCATTAACC
+
A6AAAE6EA/AAA#EAEAEEEAEEEE/EAEEEEA<E#AAAAEEEE6E<E//EE6EAEEAEAEE/AAA6E//E/EE/A<AEEEE/AEEE//<EEA/EE/E/6<<EE<EEE/<EEEEE//<AE/AAEA/

@SRR20708369.2 NS500749:581:HGV7VBGXB:1:11101:22488:1048 length=127
AGCTATCACGCCANAAGGCTTCCGAGTTAGTATAGGNGCAAAGTGAATACATCTCTGTAACCAGACTCAGATCAAGAAATGATACATTCTGGCCAGGTGTGGTGGCTCATGCCTGTAATCTCAGCAC
+
AAAAAE/EAAAAA#EEEEEEEEEEEEEEEEEEEEE6#AA/AEA/EEEEEEEEAEEEEAEA/EAE/AEAEEA6</6E///EE//6A//AEAEE<E/EE<E<EEEEEEEEEAEAEAAE///EEE/EE6/

```

:

# Answer
Cross-posted and answered on biostars: https://www.biostars.org/p/9583849/#9583850

> https://trace.ncbi.nlm.nih.gov/Traces/?view=run\_browser&acc=SRR20708369&display=metadata is a 10x single-cell sample.
>
> So FastQC plot you are showing us appears to be the R2 read from this sample which is 91 bp (split files, where as the example read you show above seems to be Illumina Index + R1 + R2). This is normal. Illumina index (8 bp), Read 1 is 28 bp (Cellbarcodes + UMI) and R2 is RNA (91 bp) = 127 bp.
> 1 votes

---
Tags: fastq, fastqc
---


# Question
Title: Remove/delete sequences by ID from multifasta
I have a fasta file like this:

```
>Id1
ATCCTT
>Id2
ATTTTCCC
>Id3
TTTCCCCAAAA
>Id4
CCCTTTAAA

```

I want to delete sequences that have the following IDs.

```
Id2
Id3

```

The IDs are in a .txt file, and the text file will be used to match and delete those sequences.

My output should be a `fasta` file like this

```
>Id1
ATCCTT
>Id4
CCCTTTAAA

```

But I want it with `awk` and/or `sed` and/or `bash` (no python or perl).
How can I do it?

# Answer
> 7 votes
Suppose you keep sequence names in `ids.txt` and sequences in `seq.fa`:

```
awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

```

# Answer
> 7 votes
Just today I wrote a script to do exactly this using `Biopython`. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script `filter_fasta_by_list_of_headers.py` is :

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

from Bio import SeqIO
import sys

ffile = SeqIO.parse(sys.argv[1], "fasta")
header_set = set(line.strip() for line in open(sys.argv[2]))

for seq_record in ffile:
    try:
        header_set.remove(seq_record.name)
    except KeyError:
        print(seq_record.format("fasta"))
        continue
if len(header_set) != 0:
    print(len(header_set),'of the headers from list were not identified in the input fasta file.', file=sys.stderr)

```

and usage :

```
filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter > filtered.fasta

```

P.S. it's quite easy to turn over the script to extract the sequences from the list (just the print line would have to move after line `header_set.remove(seq_record.name`)

I made a modification to use `in` and benchmarked it on normal insect genome assembly - 1Gbp, 300,000 scaffolds and I gave it lists of 5,000 and 100,000 scaffolds to filter. Regardless of the size of the list or algorithm used it took 29-30s. In other words, if LBYL is faster than EAFP, it's on a scale that totally does not matter.

# Answer
> 6 votes
If you want to learn how to do things with command-line tools, you can linearize the FASTA with `awk`, pipe to `grep` to filter for items of interest named in `patterns.txt`, then pipe to `tr` to delinearize:

```
$ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" > out.fa

```

This will work on multiline FASTA.

# Answer
> 3 votes
Using the `FastaToTbl` and `TblToFasta` scripts I have posted before, you can do:

```
FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa

```

# Answer
> 2 votes
Presume grep is OK if sed, awk are. The -A(n) and -B(n) flags give (n) lines post match. Presuming all your fasta sequence will be one line is a bit dangerous, but it works for your example. First get those IDs to be removed (in rmid.txt), then inverse grep against the initial fasta.

```
grep -A1 -f rmid.txt fasta.fa > rmfile.fa
grep -v -f rmfile.fa fasta.fa

```

The real answer is to use a script that defines a delimiter other than newline "\n", then parse for IDs, so better to use one of the languages you don't want to use...

# Answer
> 1 votes
This is what **`samtools faidx`** is intended for.

It needs to be called twice.

1. First, to create a FASTA index (`*.fai` file):

```
samtools faidx input.fasta

```
2. Secondly, to filter the FASTA file with the help of the index:

```
samtools faidx -o output.fasta input.fasta ids…

```

`ids…` are the IDs to retain, as individual arguments separated by space. If you’ve got a file `blocklist.txt` with IDs you want to *discard* (one per line), you first need to invert this, after having created the index (using Bash syntax):1

```
remove_ids=($(awk '{print $1}' input.fasta.fai | grep -v -f blocklist.txt))

```

… and then we can invoke the command as

```
samtools faidx -o output.fasta input.fasta "${remove_ids[@]}"

```

What’s nice about `samtools faidx` is that it also works with BGZF-compressed FASTA data — i.e. `fa.gz` files (but *only* if they were compressed using `bgzip`, which ships with Samtools, and the index must be generated on the compressed file). However, it will always write FASTA. If you want gzipped FASTA output, you need to manually stream the result into `bgzip`:

```

samtools faidx input.fasta.gz "${remove_ids[@]}" | bgzip -c >output.fasta.gz

```

---

1 **Careful**, this syntax is *generally unsafe* for reading a command output into a variable; we can use it here only because we assume that we want to generate a *space*-separated list of identifiers; the generally safe way is a lot more convoluted.

# Answer
> 1 votes
Shortest way...

```
awk '!/Id2|Id3/' RS=">" ORS=">" fasta.fa

```

`!`= Not match the pattern `Note:` Sequence doesn't need to be in a single line!

# Answer
> 0 votes
Firstly, you have to install `biopython` using the cmd `pip install biopython`.

Suppose you have two files `id.txt` and `sequence.fasta`. The solution file name is `seq_rm.py` which contains the following codes.

```
import argparse
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--id_file', required=True,
                        help='file path of removable ids')
    parser.add_argument('--input', required=True,
                        help='path of fasta sequence')
    parser.add_argument('--output', required=True,
                        help='path of selected fasta sequence')

    args = parser.parse_args()

    with open(args.id_file, mode='r', encoding='utf-8') as file:
        ids = set(file.read().splitlines())

    records = [SeqRecord(record.seq, id=record.id, description='', name='') for record in SeqIO.parse(args.input, 'fasta')
               if record.id not in ids]
    SeqIO.write(records, args.output, "fasta")

```

Run `cmd`:

```
python3 seq_rm.py --id_file id.txt --input sequence.fasta --output selected.fasta

```

It will generate an output file named `selected.fasta`.

# Answer
> 0 votes
One way using without getline:

```
awk 'FNR==NR { a[$0]; next } /^>/ { f = !(substr($0, 2) in a) }f' ids.txt seqs.fa

```

# Answer
> 0 votes
I know you do not want a Python solution but nevertheless I would recommend it like others. Here is one way to do it using Biopython:

```
from Bio import SeqIO

ids_file = 'ids.txt'
input_f  = 'input_file.fasta'
output_f = 'result.fasta'

my_dictionary = {}
with open(ids_file, 'r') as f:
    for line in f:
        my_dictionary[line[:-1]]='value'

with open(input_f, 'r') as input_handle, open(output_f, 'w') as output_handle:
    for seq_record in SeqIO.parse(input_handle, "fasta"):
        if not seq_record.id in my_dictionary:
            SeqIO.write(seq_record, output_handle, "fasta")

```

---
Tags: fasta, awk, shell, benchmarking
---


# Question
Title: Fetching Protein Sequences through UniPROT gives ConnectionError
I have 193000 protein interactions in CSV named `all_proteininteractions.csv` which have query protein name in the 2nd column and partner protein in the 3rd like this:

```
Query_ENSP,Query_Name,Partner_Name,Combined_Score
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,RAC2,0.4
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,PPP2CB,0.408
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,R3HDML,0.413
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,PLEK,0.422
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,RAB11FIP5,0.501
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,AP3B1,0.553
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,CYTH4,0.579
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,RABAC1,0.658
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,YKT6,0.672
9606.ENSP00000000233,ARF5,BET1,0.679

```

I am using this Python code to fetch sequences of each protein and store them in a CSV file:

```
import pandas as pd
import requests
import math

requests.packages.urllib3.disable_warnings()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Process each chunk
for chunk in csv_reader:
    # Fetch protein sequences for query and partner names
    query_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [fetch_protein_sequence_from_uniprot(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

```

First, it was throwing SSLError, so I made `verify=False`. But now, after running the script for 3 hours it's showing following error:

```
ConnectionError: ('Connection aborted.', RemoteDisconnected('Remote end closed connection without response'))

```

How can I remove this connection aborted error to fetch 386000 sequences?

# Answer
> 1 votes
Your script works as expected on the example you gave. Running your script results in a `protein_interactionswithseqs.csv` with lines like:

```
$ head -n2 protein_interactionswithseqs.csv
Query_ENSP,Query_Name,Partner_Name,Combined_Score,Query_Sequence,Partner_Sequence
9606.ENSP00000000233,ARF5,RAC2,0.4,MGLTVSALFSRIFGKKQMRILMVGLDAAGKTTILYKLKLGEIVTTIPTIGFNVETVEYKNICFTVWDVGGQDKIRPLWRHYFQNTQGLIFVVDSNDRERVQESADELQKMLQEDELRDAVLLVFANKQDMPNAMPVSELTDKLGLQHLRSRTWYVQATCATQGTGLYDGLDWLSHELSKR,MQAIKCVVVGDGAVGKTCLLISYTTNAFPGEYIPTVPGRPCALLQDAEAGLSSDVEGACWCQTQEGRMTGEAPAPPWPALRGPHHS

```

However, this is horribly, horribly inefficient! You are fetching the same sequence multiple times. With your example data, the protein ARF5 is actually downloaded 19 times, meaning this will take 18 times more time than needed. This is very likely the cause for the timeout error you describe. Now, it might be that your sequences are too many to be stored in RAM (but I would still try it), so one way around this is to download the sequences first, and then write the mappings. Something like this:

```
import pandas as pd
import requests
import sys
requests.packages.urllib3.disable_warnings()

def write_protein_seq_file(protein_name, protein_seq):
    file_name = "seqs/" + protein_name
    seq_file = open(file_name, "w")
    data = f"{protein_seq}\n"
    seq_file.write(data)
    seq_file.close

def read_protein_seq_file(protein_name):
    seq_file = open("seqs/" + protein_name, "r")
    data = seq_file.read()
    seq_file.close
    return data.rstrip()

# Function to fetch protein sequence from UniProt
def fetch_protein_sequence_from_uniprot(protein_name):
    #Declaring UniProt API
    uniport_api_url = f"https://rest.uniprot.org/uniprotkb/search?query=gene_exact:{protein_name}+AND+organism_id:9606&format=fasta&size=1"
    response = requests.get(uniport_api_url,verify=False)

    # Parse the response to extract sequence
    sequence = ""
    if response.ok:
        lines = response.text.split("\n")
        for line in lines:
            if not line.startswith(">"):  # Exclude header lines
                sequence += line
    return sequence

# Reading the input CSV file
input_csv = "all_protein_interactions.csv"
chunksize = 1000  # Batch size

# Create an iterator to read the CSV in chunks
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

# Initialize an empty list to store modified chunks
modified_chunks = []

# Get all sequences
for chunk in csv_reader:
    for protein in set(list(chunk.iloc[:, 1]) + list(chunk.iloc[:, 2])):
        protein_seq = fetch_protein_sequence_from_uniprot(protein)
        write_protein_seq_file(protein,protein_seq)

print("Finished fetching sequences", file=sys.stderr)
## Now, process the file again and write the new file
csv_reader = pd.read_csv(input_csv, chunksize=chunksize)

for chunk in csv_reader:

    # Fetch protein sequences for query and partner names
    query_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 1]]
    partner_sequences = [read_protein_seq_file(name) for name in chunk.iloc[:, 2]]

    # Add fetched sequences as new columns in the chunk
    chunk['Query_Sequence'] = query_sequences
    chunk['Partner_Sequence'] = partner_sequences

    # Append the modified chunk to the list
    modified_chunks.append(chunk)

# Concatenate all modified chunks into a single DataFrame
result_df = pd.concat(modified_chunks, ignore_index=True)

# Save the updated DataFrame back to CSV
result_df.to_csv("protein_interactionswithseqs.csv", index=False)

```

First, create a directory named `seqs`, then run the script above. It will process the csv file twice. Once to get all sequences mentioned in the file and save each one in a file called `seqs/$sequence_name` (where `$sequence_name` is replaced by the actual name of each protein. Then, it goes through the file a second time, but instead of the very time consuming download, it now fetches the sequence from the file stored on your computer. Once you are satisfied this works, you can delete the `seqs` directory.

Note that I don't have access to a Windows machine, so this was tested on Linux. Also, I am not a python coder, so I am sure this ins't very good python and there will be improvements to be made, but even on the toy example you gave, this speeds things up significantly:

```
$ time your_script.py

real    0m10.186s
user    0m0.998s
sys     0m0.076s

$ time my_script.py
Finished fetching sequences

real    0m3.180s
user    0m0.456s
sys     0m0.066s

```

As you can see above, even on the tiny example with only 19 data lines, your script took 10 seconds while mine only took 3. So try this approach (or, better, take my example and make it more pythonic, add some error checking, some try/except blocks, make sure your data are correct) and you might find the timeout error is no longer an issue. You should certainly find that the whole thing takes far less than three *hours*(!) to run.

First of all, that is a really, really inefficient way of doing it. Using your example data, you are getting the same protein, ARF5, 19 times. This makes your program 18 times slower than it needs to be. You should modify the script to download each protein *once* and then write them all out.

---
Tags: python, uniprot
---


# Question
Title: problem with feature counts of RNA bulk seq paired data- in Rsubread
I did manage to bild index file (from NCBI transcriptome) and perform alingment in Rsubread in R with my fq files. I did get BAM files as a result and no error. However i am having trouble with next step as featurecount function does not produce any result( 0%). I have tried builtin annotation (hg38)

i am wondering if BAM files are ok..

this is the code i have used

```
library(Rsubread)
index_path <- "D:/FASTQ_Raw/hum_transcr"
all_files <- list.files(path = "D:/FASTQ_Raw", pattern = ".fq$", full.names = TRUE)
common_prefixes <- unique(gsub("_[12]\\.fq$", "", basename(all_files)))

output_directory <- "aligned_reads/"

dir.create(output_directory, showWarnings = FALSE)

for (prefix in common_prefixes) {

  read1_file <- paste0("D:/", prefix, "_1.fq")
  read2_file <- paste0("D:/", prefix, "_2.fq")

  output_file <- paste0(output_directory, prefix, "_aligned.bam")

  align(index = index_path,
        readfile1 = read1_file,
        readfile2 = read2_file,
        output_file = output_file,
        nthreads = 4)
}

```

so the results/stats of individual BAMs appear ok still when i try featurecount function on BAM files i get zero....i was also using built in annotation index hg38 for attempt of counting but it does not work. I had issues with getting gtf file from NCBI. What do you recommend ?

# Answer
I don't think FeatureCounts is meant to work on an alignment to the transcriptome. It's not at all clear to me that the subread aligner, or any aligner, is meant to be used on a transcriptome.

If you have the transcriptome, I'd just use pseudoaligners like Kallisto or Salmon.
> -1 votes

---
Tags: r, rna-seq, bioconductor, featurecounts
---


title: converting mouse genes to human genes

text:
# Question
Title: converting mouse genes to human genes
I'm trying to convert mouse genes from PTX data to human genes in order to do a comparison with patient data to see what genes are being conserved. I'm using this file for the orthologs. **http://www.informatics.jax.org/downloads/reports/HOM\_MouseHumanSequence.rpt**

I read the table with:

```
gene_conversion <- read.table("HOM_MouseHumanSequence.rpt.txt", sep="\t", header=T, fill=T)

```

I have subsetted the mouse genes from the PTX seurat object to only include the mouse genes with:

```
mouse_only <- subset(scData, features = mouse_1)
head(rownames(mouse_only))
[1] "Xkr4"    "Gm37180" "Gm37363" "Gm19938" "Sox17"   "Gm6085"

```

Now I'm not sure how to convert the mouse genes in my mouse-only seurat object to match the human genes in the file.

Any help would be greatly appreciated.

# Answer
> 3 votes
Mice genes are usually the same as the human ones, but lower case with only the first character as upper (`Sox17` -\> `SOX17`). `toupper(rownames(mouse_only))` will do the trick.

Unfortunately, there are exceptions, and it is always better to relay on information from databases. The function bellow takes a vector of human genes, and uses `biomaRt` to return the corresponding genes in mouse.

```
# from Human to Mouse
convertHumanGeneList <- function(x){

   library("biomaRt")
   human <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
   mouse <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

   genesV2 <- getLDS(attributes = c("hgnc_symbol"), filters = "hgnc_symbol",
                 values = x , mart = human, attributesL = c("mgi_symbol"), martL = mouse, uniqueRows=T)

   humanx <- unique(genesV2[, 2])

   return(humanx)
}

```

You simply need to invert the 2 and you should be able to convert your vector `rownames(mouse_only)`

#---- Update using SOX17 as example

# Answer
> 0 votes
```
library(dplyr)
#create the function to convert mouse to human genes
homologs = read.csv("http://www.informatics.jax.org/downloads/reports/HOM_MouseHumanSequence.rpt",sep="\t")
ms2hs <- function(gene_list){
  output = c()
  for(gene in gene_list){
    class_key = (homologs %>% filter(Symbol == gene & Common.Organism.Name=="mouse, laboratory"))[['DB.Class.Key']]
    if(!identical(class_key, integer(0)) ){
      human_genes = (homologs %>% filter(DB.Class.Key == class_key & Common.Organism.Name=="human"))[,"Symbol"]
      for(human_gene in human_genes){
        output = append(output,human_gene)
      }
    }
  }
  return (output)
}

#apply the function
msgenes=c("Zfp192", "Ifitm2", "Mettl8", "Kcnh2", "Npdc1", "Pitpnc1", "Ooep", "Mettl8")
ms2hs(msgenes)

```

---
Tags: r, seurat, single-cell, conversion
---


# Question
Title: What is the best way to acquire protein isoform sequence alignment?
\[Update\] Thanks @terdon. To clarify my question:

I have a bunch of protein isoforms sequences (produced by the transcripts in the figure) and I want to align protein products (e.g., proteins produced by the first 3 protein-coding transcripts against each other). See there are 3 blocks missing in 5' region in BRCA2-204 compared against BRCA2-206.

This case might be very easy with any MSA softwares such as clustal omega, but there are complicated cases which may trick software into aligning single amino acids instead of being splicing aware. Instead of manually adjusting the alignments, I want to know whether I can directly download, or "translate", this figure into alignment results (e.g., fasta format or clustal format)

This figure should be visible with https://www.ensembl.org/Homo\_sapiens/Gene/Summary?db=core;g=ENSG00000139618;r=13:32315086-32400268 in the gene summary tab

\[Original\] Is there a way to directly download the sequence alignment result from Ensembl? If not, what will be the best way to align exons and introns without parameter-tuning

# Answer
> 1 votes
Here is a manual (or semi-manual, at best) approach. If you need to do this for very many sequences it won't be practical, but if dealing with one or two genes, it might be enough. The idea is to:

1. Get the protein coding transcripts
2. Get the proteins they code for
3. Align

In the page you gave, there is a "Show transcript table" button near the top:

If you click on that, you get the list of transcripts and their corresponding proteins:

You can then click on the "Show/hide columns" and select "Translation ID" so you can get the Ensembl protein corresponding to each transcript:

If you now deselct everything and only leave "Translation ID" selected:

You are left with a nice list of Ensembl proteins. Copy/paste those into a text file:

```
ENSP00000369497.3
ENSP00000499438.2
ENSP00000505508.1
ENSP00000439902.1
ENSP00000514856.2
ENSP00000434898.2
ENSP00000506251.1
ENSP00000499570.2
ENSP00000433168.2
ENSP00000499256.2
ENSP00000514855.1

```

Now, go to BioMart and retrieve the sequences. Select "Ensembl Genes" as the database and then "Human Genes":

Then, click on "Filters" (1) and then "Gene" (2), select "Input external references ID list \[Max 500 advised\]" (3), set the ID to "Protein stable ID(s) with version" (4) and paste the protein list (5):

Next, go to "Attributes" (1) and select "Sequences" (2):

Now, finally, click on "Results" (1) and you will get the sequences in fasta format (2):

I'm afraid they are too long to fit in this answer, but I posted them here: https://pastebin.com/NZycAmgD.

You can now paste these sequences into your favorite aligner and get the alignment you were after. Here, I used the online clustalΩ (the image only shows the first few lines):

You can find the whole alignment in clustal format (text) here.

---
Tags: sequence-alignment, proteins, ensembl, multiple-sequence-alignment
---


# Question
Title: How many false positive duplicates are marked using just the position of first unclipped base?
In the popular `picard MarkDuplicates` tool, a read is marked as a duplicate if it has the same position as another read starting from their first unclipped base in the 5'-direction. The implication is that if two reads have this property, it must be much more likely they arose from some "artificial" duplicates (PCR, optical, etc.) rather than due to sheer coincidence that two independent fragments of DNA happened to have this property and were sequenced in the same run.

Does anyone have links to studies or other data that attempts to measure how often this mechanism for labeling duplicates fails? I.e. how frequently does a read get marked a duplicate in this way "incorrectly," meaning it was actually just a coincidence and really did correspond to an independent fragment that should be considered in your analysis.

I imagine this answer can vary quite a bit based on the sequencing conditions: library prep, machine, etc. I was hoping for at least a ballpark estimate in the case of some "standard" modern short read sequencing methods. But it should be possible to study this in particular cases, e.g. using UMIs, doing PCR vs non-PCR runs, etc.

# Answer
> 6 votes
There's a tool in Picard (full disclaimer: I wrote it...) called CollectIndependentReplicateMetrics which performs the test you are looking for, I think.

It measures the percent of duplication that is "independent", i.e. which arose from randomly sampled original fragments that *happen* to be duplicates, rather than due to technical (optical, PCR, bridge-amp...) reasons.

The way it works is by looking at duplicate-sets of size 2 (and separately, 3) over heterozygous sites in the bam (which is assumed to be a diploid sample) it classifies the sets into "homogeneous" and "heterogeneous" based on the alleles of the hetrozygous variant. From the relative sizes of these two classes, one can work out the percent of duplication that is "independent" since independent duplication will give rise to heterogeneous sets (sometimes), while dependent duplication will *never* do that.

The assumptions are quite simple actually:

* a technical duplicate will have the same allele
* an independent replicate will have the same allele 50% of the time

So the rate of independence is 2x the rate of heterozygosity, *within* duplicate sets of size 2.

For sets of size 3 the math is a little more complex, so I wrote it up and docuemnted it.

Here's the verbatim explanation I wrote at the time (in the source-code of the *Metric* file...):

```
Explanation of calculation of independent replication rate from heterozygosity rate (within triplicate sets):

We assume the there are two types of replication events:

- those that are "independent", such that we just happen to get 2 fragments from the exact same region
(These get a random allele so effectively change the allele with probability 0.5), and
- those that are "artifactual" (do not change the allele)

We skip sets that have unexpected alleles as they do not fit our model. In the following we use the term "duplicates" to indicate that
the read-pairs would be marked as duplicates, not that they actually are technical duplicates.

To reach a triplicate set, 2 replication events are required so there are the following options, assuming
that the independent replication rate is x (thus the artifactual is 1-x). We assume a diploid organism with no bias towards either allele
in heterozygous sites, so an idependent replication will result in the other alleles in half the cases:

0 -> 0, 1 happens with probability x/2, therefore
0 -> 0, 0 happens with probability 1-x/2

(This is the explanation that is required for calculating the independent replication rate from doubleton sets...quite simpler)

Without loss of generality we assume that we "start" with allele 0 and that 1 is the other allele.

Each of the resulting alleles ("First" or "second") can replicate (each with probability 0.5) so we get:

from first row:
0=>1 1 (0 replicate to a 1) with probability  x/2 * 0.5 * x/2
0=>0 1 (0 replicate to a 0) with probability  x/2 * 0.5 * (1-x/2)
=====subtotal = x/4

0 1=>0 (1 replicate to a 0) with probability x/2 * 0.5 * x/2
0 1=>1 (1 replicate to a 1) with probability x/2 * 0.5 * (1-x/2)
=====subtotal = x/4

from second row:
0=>1 0 (first 0 replicate to a 1) with probability (1-x/2) * 0.5 * x/2
0=>0 0 (first 0 replicate to a 0) with probability (1-x/2) * 0.5 * (1-x/2)   <======= Homogeneous set
=====subtotal = (1-x/2)/2

0 0=>1 (second 0 replicate to a 1) with probability (1-x/2) * 0.5 * x/2
0 0=>0 (second 0 replicate to a 0) with probability (1-x/2) * 0.5 * (1-x/2)  <======= Homogeneous set
=====subtotal = (1-x/2)/2

total is x/2 (from first two sub-totals) + (1-x/2) (from last two sub-totals) = 1

We differentiate between a heterogeneous result (with 0's and 1's in the set) and homogeneous results
(all zeros, since we assumed WLOG that we start with 0)
The probability of a homogeneous set is therefore the sum of the two only homogeneous results

P(hom) = (1-x/2) * 0.5 * (1-x/2) + (1-x/2) * 0.5 * (1-x/2) = (1-x/2)^2 = y

where y = P(hom) is the rate of homogeneity within triplicate sets.

Solving for x we find that 2 * ( 1 - sqrt(y) ) = x.

```

---
Tags: sequence-alignment, reads, markduplicate
---


# Question
Title: conducting a ttest on a board of 1000 genes in 50 peope divided into to groups(A and B). pvalues are different each time
I am struggling in conducting at test on a board of 1000 genes in 50 individuals divided into to groups (A and B). I have a dataset of 1000 genes of two groups of people (total individuals=50) group A = 30 (healthy) individuals & group B = 20 sick) individuals

I created a matrix where rows are genes and ind are columns. now I need to perform a t-test for all the genes and find the p values of my alternative hypothesis (not the null one): i)A\B and A=!B(not equal). threshold is 0.05 after FDR (which I must include)

I am struggling with the code since t-test differ on sensitivity and each of my approach returns a different result. The dummy data is as follows:

```
import pandas as pd
import numpy as np

# Generate data for group A (healthy individuals) & group B (sick individuals)
A = np.random.rand(30, 1000)
B = np.random.rand(20, 1000)
# Create dataframe for group A & B with genes as columns
Adf = pd.DataFrame(A, columns=[f'Gene_{i}' for i in range(1, num_genes + 1)])
Bdf = pd.DataFrame(B, columns=[f'Gene_{i}' for i in range(1, num_genes + 1)])

# Add a column indicating the group
Adf['Group'] = 'A'
Bdf['Group'] = 'B'

df = pd.concat([Adf, Bdf], ignore_index=True)
df.head()

```

Any ideas?

# Answer
Could it be that you are getting different results because you are not setting a seed on your random number generator?

The recommended way to make sure you get the same results when using numpy to generate dummy data is:

import numpy as np
# make random number generator with seed
rng = np.random.default_rng(seed=42)
rng.random([4,4])
#Out:
#array([[0.77395605, 0.43887844, 0.85859792, 0.69736803],
#       [0.09417735, 0.97562235, 0.7611397 , 0.78606431],
#       [0.12811363, 0.45038594, 0.37079802, 0.92676499],
#       [0.64386512, 0.82276161, 0.4434142 , 0.22723872]])

```

> 1 votes

---
Tags: python, gene, p-values, matrix, t-test
---


# Question
Title: How to calculate pharmacophore features?
I asked this question already in the matter modeling forum (here). But i hope this is a better place because the topic is more biological.

I want to use the enzyme ß-Glucuronidase with the pharmit database. I uploaded my file for the receptor from pdb file. How can I get the pharmacophore features? So one thing is the enzyme is quiet popular, so is there any database with the features? I tried it with a publication but the given link is not working. It checked the site from the author but I did not found the database there.

And if there is no database or maybe there is one but the enzyme is not included, how can I get the pharmacophore features for my calculations?

# Answer
> 1 votes
Just put a know ligand inside the active site of your pdb , the webserver will calculate the Features of the ligand/pharmacophore then you can select the ones you want and the database where to search for molecules with the same Features. All the software used to crate a Features file seems to be not for free or really old (not sure if you can trust that code anymore) if you need some particular Features , you'll have to draw a FAKE LIGAND by yourself that carries them ... but thats its just my naive idea – –

se my resul using E.coli 3LPF pdb carrying **Z77**: 1-\[(6,7-dimethyl-2-oxo-1,2-dihydroquinolin-3-yl)methyl\]-1-(2-hydroxyethyl)-3-(3-methoxyphenyl)thiourea ligand and selecting first 4 Pharmacophire Features extracted by the webserver from it and running them against the CHEMBL32 database, first hit selected.

**EDITED**

The instructions above didnt work for me , loading `3lpf.pdb` from my PC to the webserver for some reasons doesn't get the Z77 ligand. So the only option is to draw the Features by yourself using the Features menu, played a little bit with it but not well enough to get results , also the tricky part here would be to set the right coordinates for your ligand-Features, but overall the Tool looks very good, really fast and with a nice intuitive interface.

Still would like to get ligands from PC uploaded pdb file too , so that I could test my Fake-Ligand theory.

**RE-EDITED** :

found a third way :

load your pdb from your pc , save the session , modify the `pharmit.json` file you save adding the FAKE-LIGAND. To have an idea of how to put your FAKE-LIGAND inside the `json` file , have a look at a saved session from a session started suppliying a pdb code instead of a file. In my experiment I loaded the modified Z77 ligand stiched to a session started from a downloaded `3lpf.cif` and saved to a `.json`

auchhh !!!!! it does load the ligand but doesn't create the Features , don't know if it is because the ligand is cut in 2 pieces or because it just doesnit work like when I uploaded a `.pdb` containing a ligand.

At least you have a way to position the Features you'll add by hands afterwards.

\***think that this one wont never works because the Features for the ligand are saved in the .json (first entry in the dictionary) the webserver calculates them once and then save them to be reloaded); yes it does do that , I can delete the ligand entry from .json and reload the session getting the Features but not the ligaand.

**FOURTH EDIT (last one I swear)** :

Tried new approach: loading Features menu uploading Z77 from https://www.rcsb.org/ligand/Z77 page, both model and ideal molecule as `.sdf` file. The ligand gets uploaded and Features calculated by the webserver , obviously the model molecule doesn't get located in the active site of the molecule, so I guess you could draw yor FAKE ligand using something like **PyMOL** exactly where you want it to be and then upload it with the load Features menu.

My uploaded `Z77_ideal.sdf` pic:

So that would provide ground for my FAKE-Ligand approach

---
Tags: docking
---


# Question
Title: What is a common number of cells per cell type in single cell RNA Seq?
Good afternoon,

I just started working with Single cell RNA Seq and I am trying to understand what is an average / possible number of cells per sample.

I am looking at a dataset of 7 patients with kidney tumour samples (please find attached). Would these number of cells be possible? I basically want to double check whether I didnt do any mistakes in my pre-processing.

Thank you

# Answer
It depends on many factors, both biological and technical, and it's impossible to answer in general.

The starting question is which technology was used, and how many cells were used in input. For example, in a SMART-Seq experiment with one cell per well (in a multi-well plate), you should know how many cells you have (one per well). For a 10x experiment, the person running the experiment should have measured the cell concentration in the suspension before loading the chip, and targeted a particular number of cells (typically 10,000). In a combinatorial barcoding approach (e.g. Parse), you'd also have a measurement beforehand so as to load a particular number of cells per well. So, you should already have an idea of the total number of cells.

Then, the first step of the analysis is typically some QC. In barcode-based approaches (including 10x), you'd make a knee-plot to estimate the actual number of sequenced cells in your sample.

In the clustering and annotation steps, you'd start distinguishing and annotating groups of cells, which hopefully do correspond to cell types. The algorithms are not perfect, so the clusters don't necessarily match cell types perfectly. For example, a rare cell type will easily be lumped together with another, similar, cell type. For a well-represented cell type, the algorithms have more chance to capture subtypes or cell states.

And of course there is the biology, and your knowledge of it. How many astrocytes are there in the kidney? I assume any number above 0 is worrisome. I'd expect the cell types in the data to roughly reflect the cell types in the tissue (but not exactly, there can be all kind of bias).

There is a question of depth: if you have few cells in a cluster, and few reads per cell, the gene expression data in that cluster is not very reliable. The actual number depends on the technology and read depth. I would say for 10x I would prefer at least 100 cells, and would not trust a cluster with \<10 cells.

So for your question, you need to use your knowledge of the kidney and of the experimental design to decide. You can also compare to published scRNA-Seq experiments of the same tissue. You can also compare between patients to see if there is something troubling (for example, why does patient 1 have epithelial cells but no adipocytes, and opposite for patient 7? It might not be a problem, but only with context can you tell). And importantly, you do need to check that the annotation does match known gene expression patterns.
> 0 votes

---
Tags: scrnaseq
---


title: How to adjust cutoffs in Mutect2

text:
# Question
Title: How to adjust cutoffs in Mutect2
I'm currently performing some initial analyses on targeted Illumina sequencing data generated with DNA from diseased tissue. My hypothesis is that one particular gene might contain somatic mutations in this disease. I'm currently running Mutect2 in tumor-only mode, using default settings.

As I have only found very few mutations (this result might well be plausible), I'd like to adjust the Mutect2 cutoffs to increase the sensitivity and see whether the result changes. For instance, I'd like to reduce VAF, depth of coverage and QUAL. Unfortunately, I can't work out which of the many arguments have to be modified and how: https://gatk.broadinstitute.org/hc/en-us/articles/360037593851-Mutect2

Can you help me out? Thanks!

```
gatk Mutect2 \
-R $fastaFile \
-I $inputFile \
--germline-resource $germlineResource \
-O $resultDir/$sampleName.Mutect.vcf.gz

```

# Answer
Mutect2 uses a statistical model based on the available reads. Unlike with some other tools, there is not threshold based on depth or VAF or QUAL. The main statistic that *does* govern emission is the LOD (Log ODds) which is the log10 of the odds ("ratio") that the variation being inspected is a tumor relative to the probability that it is a technical artifact. The default threshold for this is 3.0 (corresponding to a odds ratio of 1:1000).

Perhaps you'd like to see lower probability results, in which case you should lower the threshold:

Using `--tumor-lod-to-emit 2.0` should give you many more results, many of them might be wrong. The LOD isn't well calibrated, so I'd take it with a grain of salt....**validate** the results you get.
> 4 votes

---
Tags: gatk, somatic, mutect2
---


# Question
Title: What does a question mark ("?") mean in Picard metrics files when I expect a number (integer, float, etc.)
Does anyone know what conditions must be true for Picard to put a `?` where I expect a number?

Here is an example from the output of the tool `CollectWgsMetrics`:

```
❯ cat sample.wgs_metrics.txt | head -n8 | tail -n2 | verticalize | grep FOLD
FOLD_80_BASE_PENALTY    1.33476
FOLD_90_BASE_PENALTY    2.288159
FOLD_95_BASE_PENALTY    ?

```

Code searching (link) didn't reveal anything, and googling isn't returning anything obvious either.

---

*\[Verticalize is a little CLI tool I've forever been using by Pierre Lindenbaum: `verticalize`. Easy to compile and I toss it in my `~/bin` every time I get a new computer.\]*

# Answer
> 5 votes
Picard (actually HtsJdk) uses a formatting class (via the MetricsFile class) that emits "?" for NaNs (and Infs).

In your case, there seems to be a division-by-zero error during the calculation of the metric (likely due to zero coverage at some small percentage of bases in your input). This resulting Nan is emitted as a "?".

---
Tags: picard
---


# Question
Title: Variant allele frequency versus population frequency
In gnomAD related article

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC9160216/#:~:text=Allele%20frequency%20is%20calculated%20by,can%20differ%20substantially%20between%20positions

I road this section

> Allele frequency is calculated by dividing allele count by allele number, hence, allele frequency represents the frequency of confidently sequenced haplotypes that carry the allele in question (because coverage and sequencing quality varies across the genome, the allele number can differ substantially between positions). Allele frequency is not equivalent to the percentage of individuals that carry the allele, but is a suitable value for expressing the frequency of a variant in the general population. The number of individuals carrying a variant will depend on the number of heterozygous and homozygous individuals but can be calculated from the data provided in the variant table.

Does this mean, for getting the frequency of each allele in gnomAD, one could use `Allele frequency (AF)` column or must calculate population frequency?

# Answer
gnomAD provides a global Allele Frequency as well as AF values for most subsets. gnomAD does not provide individual genotypes (as far as I know), so calculating custom AF for non-precomputed subsets might not be possible.

Since you're interested in the NFE (non-Finnish European) sub-population, the `nfe_*_AF` fields might be of interest to you.

---

As far as calculating the fraction of non-WT samples from `num_homozygous_samples`, `AC`, and `AN` goes, your formula seems spot on. Here's my deconstruction of the formula:

Your formula is: `(AC - num_homozygous_samples)/(AN / 2))`

`AN/2` is the number of samples (assuming all samples are diploid). `AC - num_homo` gives you `2 * num_hets`. Dividing this by number of samples gives you frequency of samples with at least one ALT allele.

For example, let's take a dataset of 200 diploid samples where 100 are WT, 75 are het and 25 are hom-alt for a biallelic variant.

Since there are 200 diploid samples, `AN` is 400.

`AC` for the ALT allele is `(75 * 1) + (25 * 2) = 125`.

Since 100 are WT, it is known that fraction of samples with mutation = `1 - (100/200) = 0.5`. This number can be arrived using your formula:

```
(AC - num_homozygous_samples) / (AN/2)
 = ((125 - 25) / (400/2))
 = 100 / 200
 = 0.5

```
> 3 votes

---
Tags: allele-frequency, gnomad
---


# Question
Title: Converting .tsv to .meg for visualization
This is an add-on to the previous question I asked here.

I was able to deploy a "chunking" strategy to my kmer matrix and was able to get a file with the distance measurements! From what I could tell by using the `head` command was that the data looks like what it did before when I was using this approach on a smaller dataset so now I'm going to use this data to generate a tree using MEGAX.

Normally when I'd convert the `tsv` file to a `meg` file I'd manually add in the meg headers/requirements for conversion. However, I can't open the distance `tsv` file (its ~23 GB) to make these changes.

Is there an easy solution to this problem? How do you convert a `tsv` file to a `meg` file without opening the file?

EDIT:

For more context here is the format for each file.

Currently, my `tsv` file looks like this:

```
0 1 2 3 4 5 6 7 8 9 10 ... 178
1 0 # # # # # # # # #  ... #
2 # 0 # # # # # # # #  ... #
3 # # 0 # # # # # # #  ... #
4 # # # 0 # # # # # #  ... #
5 # # # # 0 # # # # #  ... #
6 # # # # # 0 # # # #  ... #
7 # # # # # # 0 # # #  ... #
8 # # # # # # # 0 # #  ... #
9 # # # # # # # # 0 #  ... #
10 # # # # # # # # # 0 ... #
...
178 # # # # # # # # # #  ... 0

```

This is a jaccard distance matrix of 178 samples. The `#` represent the distance values for each sample compared to the other. I was able to calcualte the distance using pandas and scipy: `matrix = pd.DataFrame(squareform(pdist(data, 'jaccard')))`

And I would like it to look like this:

```
#mega
TITLE: Jaccard distance of 178 samples

sample_1
sample_2
sample_3
...
sample_178

distance_1
distance_1 distance_2
distance_1 distance_2 distance_3
distance_1 distance_2 distance_3 distance_4
...
distance_1 distance_2 distance_3 distance_4...distance_178

```

# Answer
> 1 votes
Note, its a while since I've used `scipy` and I don't fancy loading an deleting for this test ... so in the absence of `pdist()` the solution would be ...

```
import numpy as np
import pandas as pd
from pathlib import Path

data = [[1,2,3],[0,5,6],[7,0,9],[10,11,0]]
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan'
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file = Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

```

**output**

```
#MEGA
a
b
c
d

0
1.0 0.0
5.0 7.0 0.0
9.0 10.0 11.0 0.0

```

That works. Keep in mind where it says `arr.reshape(3,-1)` the `3` should be the number of columns -1 in your matrix: **in your case this appears to be 177 (178 - 1)**.

The gist is to use `np.tril`, which doesn't exists in `pandas`. Hence don't use `pandas`.

At a guess with `scipy` then it would be ...

```
import numpy as np
import pandas as pd
from pathlib import Path
from scipy import pdist

tmp = [[1,2,3],[4,5,6],[7,8,9],[10,11,12]]

data = pdist(tmp, 'jaccard')
arr = np.tril([data]).astype('float')
arr[arr == 0] = 'nan'
arr = np.where(np.isnan(arr), '', arr)
columns = ['a','b','c', 'd']
file =Path('/Volumes/Data2/test/test1.meg')
if file.is_file():
    file.unlink()
with open(file, 'a') as fout:
    fout.write("#MEGA\n")
    for col in columns:
        fout.write(col + "\n")
    fout.write("\n0\n")
    np.savetxt(fout, arr.reshape(3,-1), fmt='%s')

```

You might need a quick `perl -p -i -e 's/[" "]+/ /g' test1.meg` or even `perl -p -i -e 's/[" "]+/\t/g' test1.meg` to set the file for MEGA. I dunno why but `\s` or `" "` in the transform part of perl pie isn't working for me.

---
Tags: format-conversion, k-mer, megax
---


title: Low custom tdtomato gene content

text:
# Question
Title: Low custom tdtomato gene content
I have a set of scRNA-seq samples expressing TdTomato, which has high content in microscope. I followed the 10x cellranger pipelines to finsh the work, my procedures are as follows:

1. added TdTomato on an additional chromosome:
2. added an entry in the gtf
3. cellranger mkref
4. cellranger count

I got the result, which includes:

But strangely, the tdtomato expression content in my three sample are all very low! And cells with tdtomato are also very low and disperse!! Take sample 1 as example, you can see the picture:

It's not normal, but I can't find the reason. Following the 10x genomics help documentation noticed that by executing the CellRanger filtering steps, maybe my tdtomato reads are filtered too much, but because I'm not very skillful at programming and I'm first time to do RNA-seq (and scRNA-seq), I just don't know how to analysis the bam file as suggested by the official website and find the real reason of the problem:

If you can show me the code need to do, I will be very grateful!

# Answer
> 3 votes
A few suggestions:

1. Manually inspecting the BAM file with a genome browser like IGV, as @GWW suggested and checking the qualities of the alignments in your BAM/SAM file (SAM is the human readable version of BAM, there are tools to convert one to another);
2. Using another single cell analysis tool like Seurat so that you can manually set thresholds for QC and not lose to many (tdtomato) expressing cells at the QC step (I don't know about your experiment but if transfected, cells could be at a worse condition than non-transfected ones and hence might not pass QC);
3. Double checking your tdtomato reference sequence, I don't know about tdtomato but for example GFP has flavors like EGFP, similar but different sequences, might explain low alignment score if you observe in step 1;
4. Taking the limitations of scRNA-seq (especially droplet-based) technologies into account. The figure below is from the CITE-seq paper (from supp info), *Simultaneous epitope and transcriptome measurement in single cells* (Stoeckius et al., 2017), showing that the protein - RNA measurements of a given gene might show poor correlation.

# Answer
> 2 votes
The solution to your problem probably is adding the full mRNA sequence of your transgene to the reference (as also suggested by acrux).
I had a very similar problem recently when tdTomato expression was only detectable in a couple of cells. The transgene was introduced into the cells using a lentiviral vector and therefore the mRNA had the tdTomato sequence followed by a WPRE element. In the picture below you can see that a large number of reads is mapping to the WPRE element in the 3' UTR (blue region) and and almost none to the tdTomato CDS (red region).

In general, the untranslated regions are considered exons because they are part of the final mRNA. So if you have the fasta sequence beginning at the transcription start and ending at the polyadenylation signal, this is the full mRNA. Then you just need to update the GTF with the longer full mRNA length and you are all set.

# Answer
> 1 votes
We had the same issue. Turned out, our tdTomato transgene was followed up by WPRE and BGH polyA signals. These sequences are quite large, and their presence precludes sequencing into the tdTomato coding sequence when using the 10x genomics platform. Incorporating these into our build, allowed us to detect tdTomato+ cells.

# Answer
> 1 votes
I am unable to respond to other's comments as I never usually comment, so this is a new account. However, I spent the last couple days really diving into the BAM alignments when mapping scRNAseq reads to the human genome, but altering the GTF and FASTA using different annotations of my CAR construct (some very liberal allowing upstream and to the poly A tail sequence, and some less liberal)

MY MAIN CONCLUSION!:

DO NOT ADD THE BGH POLY A SIGNAL TO THE GENE ANNOTATION! it is highly homologous to the human beta globin subunit, and I am finding almost none of the reads actually mapping both to the adjacent CAR construct and the BGH POLY A subunit.

I am still figuring out how to handle the CD3zeta, CD28, and the OX40 in the CAR construct being multi-mapped (I will test just using a small portion of the CD28 transmembrane domain and then using just the chimeric antigen portion for mapping)

BUT, BIGEST TAKEAWAY, if you are going to attempt to add the poly A tail or any other genomic annotation not strictly relating to your gene of interest look at the bam and ensure that the reads mapped actually bridges into your gene. (as I also had residual genomic DNA from human beta globin upstream of my CAR construct that showed splicing over the CAR construct and aligning downstream of the CAR to the POLY A tail sequence)

-for viewing the BAM i am using IGVtools, and for the creation of the GTF and the FASTAs that I concatenated to the human reference genome I used tears \[not a software package\] ;)

# Answer
> 0 votes
you might want to try adding the 3'UTR sequence of you tdtomato construct used in your cells. that should do the trick.

---
Tags: scrnaseq, single-cell, 10x-genomics, cellranger, tdtomato
---


# Question
Title: Dog (Canis lupus familiaris) MT gene annotations
I have an RNA-Seq and scRNA-Seq dataset from dog samples and pseudo-aligned them to the dog reference transcriptome (Ensembl release 109). When filtering for mitochondrial gene percentage in scRNA-Seq I couldn't find the "MT" annotation in this new release. Does anyone have any information about this? Did they not include it or am I doing something wrong somewhere (the MT annotations are present in the CanFam3.1 genome built but not in the newer ones). I could not find anything on the Ensembl release details or blog.

Thank you :)

This is a part of the code that I used to search for the mitochondrial genes:

```
searchFilterOptions(mart=ensembl_canis, filter="chromosome_name", pattern="MT")

```

# Answer
Currently, there is no available mitochondrial sequence for the ROS\_Cfam\_1.0 assembly: https://www.ncbi.nlm.nih.gov/datasets/genome/GCF\_014441545.1/
> 3 votes

---
Tags: sequence-annotation, ensembl
---