title: Sequence alignment using BWT

text:
# Question
Title: Sequence alignment using BWT
**My Problem:**

Skipping some specific background, what I want to do is **judging whether some soft-clipping sequences are the same**, which may result by the same SV event. Colored bases in Fig.1 is an example of soft-clipping sequences.

I use `BWA -MEM` as the aligner. And I am trying to develop an SV caller actually, what I want to do here is refining PE(Paired-Ends) signal. If it is SR(Split-Reads), I don't need to do this, because of SR and its supplemental alignments coming from the same read. But PE signal generated from some cluster procedures need some filters. And I think K-mer is not a good way to do this which adapted by some SV callers. For me, as a bioinformatic fool, imitation and implementation may be a good learning process. It seems still a long way to go.

**My solution**:

I suppose this is a sequence alignment problem. Then I use FM-indexing based on BWT to do this. My inexact/approximate matching implementation is referred from BWA.

**My workflow**:

1. Extract these soft-clipping sequences
2. Choose the longest one as the reference
3. Iterate others to match with the reference, if the number of unmatched is large than 1, it means these sequences are not the same.

---

* Question 1: Is that sounds reasonable to you? Any suggestions?
* Question 2: If it is OK, how should I implement **the penalty scoring scheme**?
Here, because these sequences should come from the same SV event, the main differences(Mismatch, InDel) should arise from sequencing errors, not small genomic variants. So I want to set a higher gap open penalty and even higher gap extend penalty, like

mismatch\_penalty = 4 gapopen\_penalty = 11 gapextend\_penalty = 20 alignment penalty cutoff = 30

The later means that one alignment with two more gaps will be considered as spurious alignment.

How do I optimize these parameters?

# Answer
You just need to check the flag for suplementary mapping. If there is suplementary mapping then use the coordinates; if they are all the same the seq is the same and you just have one event, if different reads have different suplementary mapping locations then there's multiple events. The problem is more complex than this though so just use an existing tool.
> 1 votes

---
Tags: sequence-alignment, bwa
---


title: Building small molecules de novo

text:
# Question
Title: Building small molecules de novo
Can anyone recommend a protocol for building the 3D structure of a small molecule *de novo*, given only a 2D rendering? So far my idea has been to build the molecule using the fragment builder in PyMol, then make conformers using Frog2. But I have no idea how realistic the starting structure is from PyMol and how exhaustively the conformers are sampled in Frog2.

The downstream application of this would be to dock the ligand into a protein binding pocket. If anyone has a more rigorous protocol for making sure all the conformational space is fully sampled, it would be greatly appreciated!

# Answer
Most universities have a ChemDraw site lincence, so you may want to check that for drawing a non–3D-embedded molecule. Alternatively there are many other programs. PyMOL is not really a good choice. Also there are online drawing tools everywhere, where you can download a SMILES (see below).

Chemical formula is the wrong way to go, unless you are "enumerating". The chemical formula of glucose is `C6H12O6`. This is the same for all hexoses. What you want is a SMILES. This string gives the connectivity of the atoms. For D-glucose it's `OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O`, the `@` symbols are chiral centres. In Cheminformatics making a list of all possible compound with n-carbons is a big deal, this is called enumeration and you likely don't have to reinvent wheel. Furthermore, there are datasets that are filtered by the Lipinski rule.
> 1 votes

---
Tags: docking, computational-biochemistry, small-molecules
---


# Question
Title: Splitting a GenBank file into smaller files
I have a genbank file that I have been trying to run in HMMER using the following command:

```
hmmscan --domtblout NC_000962.3.final.gbk.out -E 0.000001 prot_profile.hmm NC_000962.3.final.gbk

```

However, since the file is too big for HMMER, I get the following error:

```
Fatal exception (source file p7_pipeline.c, line 697):
Target sequence length > 100K, over comparison pipeline limit.
(Did you mean to use nhmmer/nhmmscan?)
Aborted (core dumped)

```

I need to split this file into smaller files for the HMMER run. The problem is, I would need the *ORIGIN* field (containing the nucleotide sequence) at the bottom of the file split into relevant chunks into the subfile as well. This is where I am getting stuck.

I have already looked at *GenBank Parser* and *seqretsplit* to no avail.

I would greatly appreciate some help with this. Thank you.

# Answer
You can use NCBI EntrezDirect to download data in GenBank format for a specific region of the sequence as follows:

```
efetch -db nuccore -id NC_000962.3 -format gb -seq_start 1 -seq_stop 99000 > seq.gb

```
> 1 votes

# Answer
I don't have time to give you a polished fully functional version but this will get you started. The next step will be a bit harder as you will need to capture the start of the first CDS in the chunk and end of the last CDS instead of just hardcoding as I have here.

```
fout = open('/path/first3.gbk','w')
for record in SeqIO.parse('/path/NC_000962.3.final.gbk','genbank'):
    for i, feature in enumerate(record.features):
        new_record = record
        new_record.seq = new_record.seq[0:4437]
        new_record.features = new_record.features[:3]
        if i > 2:
            SeqIO.write(new_record, fout, 'genbank')
            fout.close()
            break

```
> 1 votes

---
Tags: genbank, hmmer
---


title: Clustering RNAseq fold-change data

text:
# Question
Title: Clustering RNAseq fold-change data
I have a dataset that looks something like this

|  | Treatment1 | Treatment 2 | Control |
| --- | --- | --- | --- |
| Sample1 | 3.23. | 0.87. | 1 |
| Sample2 | 1.71. | 1. | 1 |
| Sample3 | 2.88. | 5.65 | 1 |
| Sample4 | 0.77. | 1.34 | 1 |

The numbers describe a fold change with respect to the control (where the control is always set to 1).

Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about *just* fold change?

By extension, does it also make sense to compute something like Euclidean distances (e.g. fold change vector of Sample1 vs. Sample2) for clustering fold change data? I think clustering on counts/TPM makes more sense, but wasn't sure about doing it on fold changes.

Something tells me that it's not possible because the fold changes are relative values as opposed to measured numbers. I've opted for looking at correlations (e.g. Spearman) for now though; are there alternative metrics/strategies to consider?

If there are some good papers talking about clustering and manipulating fold change data I'd be equally grateful, thanks!

# Answer
> 1 votes
> Does it make sense to compute values such as averages, sums, and differences (e.g. the average fold change for treatment1)? And yes, wary that doing things on log fold-change is not a thing, but how about just fold change?

This is problematic, if you want to compare treatments then compute the actual fold change and its significance. You can't really base anything off of the fold-change (log or otherwise) alone.

> does it also make sense to compute something like Euclidean distances for clustering fold change data

It's not unheard of to do this. The problem is you lose all information about how reliable those fold-change estimates are. This is something maintained when you cluster sample-level information (e.g., TPM or robust log transformed data). Not that using variance-stabilized data for clustering will work better, since it tamps down on how much highly-expressed genes influence your clustering (assuming you don't want that to occur).

Have a look at the DESeq2 vignette for some nice examples of what good practices are.

---
Tags: rna-seq, differential-expression, clustering, fold-change
---


# Question
Title: Is it possible to filter contaminated reads for raw PacBio sequences (not HiFi reads) before assembly?
De novo genome assembly for non-model organisms face the issue of bacterial contamination. For assembled contigs with mostly bacterial-like sequences (based on BLAST search), the entire contig can be dropped off. But for a long contig with small fragments of bacterial-like sequences, one way to solve it is to break the contigs into smaller pieces and drop the bacterial-like pieces. It would work better if the raw reads could be removed BEFORE assembly. This was commonly done for Illumina reads, but for PacBio long reads (non-HiFi reads) with a higher error rate, is it possible to filter contaminated bacterial sequences based on sequence similarity, like BLAST, Kraken? Of course, the potentially contaminated region should be long enough to drop the entire sequence, for PacBio long reads. Maybe 30%? Any suggestion of the processing details?

# Answer
> 1 votes
I don't think that in principle the issue is any different for PacBio and Illumina; it's just that the process might be noisier with higher error rates. Note that PacBio errors aren't randomly distributed substitutions and probably many (most?) short-ish k-mers (\<30nt) will likely be intact. I would suggest doing k-mer filtering of your reads as suggested here or here or here.

Alternately, use your favorite alignment-based technique, if you have a good idea of what your contaminant is. In principle PB reads will be *better* filtered by this method than Illumina as local alignment will be able to figure out the errors quite well with longer queries.

If you go the route of trying to do alignment to some set of contaminant references, the obvious answer would be to use minimap2 as your aligner as this is exactly the task that it is designed for. You can then use samtools to filter out anything that maps well.

Here is a related question/answer that also addresses the contaminant issue and suggests using Blobology, which I'm not familiar with but may be helpful (granted it also seems to be Illumina-centric). Genome-size estimation using k-mers on PB reads appears to be pretty accurate too, so I think that it's likely ok to work with raw PB k-mers.

# Answer
> 0 votes
You can use the kraken database. It will classify your reads based on what organism it matches with. Then can write a custom script to extract the reads you are interested in for downstream analysis.

---
Tags: genome, assembly, long-reads, pacbio
---


# Question
Title: Looking for a tool to find 16S rRNA in hundreds of genomes
I have hundreds of bacterial genomes and I would like to know which is the best informatic tool in which I can run all genomes on one go and get a 16S rRNA gene sequence of each genome as output.

Could anyone give me some advice?

# Answer
You can use Barrnap

```
barrnap -k bac --threads "$(nproc)" \
        -o output_rrna.fna < input_file.fna \
        > output_rrna.gff 2> /dev/null

```

Put the file in one line of sequence per fasta id and get sequence with 16S annotation.

```
perl -pe 'if(/\>/){s/\n/\t/}; s/\n//; s/\>/\n\>/' output_rrna.fa | \
    perl -pe 's/\t/\n/' | tail -n+2 | cut -d\: -f1 | \
    grep -A1 16S > output_16S.fasta

```

Check the `output_rrna.gff` and make a BLAST of the sequence to be sure about the predicted 16S.
> 2 votes

# Answer
You can perform genome annotation using tools like prokka or PGAP. Using some python scripts you can extract the predicted 16S and blast to confirm.
> 0 votes

---
Tags: genome, bacteria
---


# Question
Title: Getting Unique Identifier List for GEO Datasets NCBI
**AIM:** Download "Unique Identifier List" for the following query from GEO DataSets.

**Query:** *("Expression profiling by high throughput sequencing"\[DataSet Type\] AND ("Homo sapiens"\[Organism\] OR "Mus musculus"\[Organism\] OR "rattus norvegicus"\[Organism\])) AND ("2020/01/01"\[PDAT\] : "3000"\[PDAT\])*

which means, all RNASeq studies deposited on GEO in the year 2020 for humans, mice or rat.

Problem: I need the GSE ID list for ~9k datasets, but while trying to download the list of ids, it loads to a blank page and nothing happens. Also, clicking on "Next Page" gives error. I have been trying for the last 3-4 days but it doesn't work.

**Steps to generate file:**

*"Send To" -\> "File" -\> "Format" (Unique Identifier List) -\> "Sort By" (Default Order) -\> "Create File"*

# Answer
> 1 votes
You can use Entrez Direct for this. The following returns Unique Identifiers which are just bare integers.

```
$ geo_query='"Expression profiling by high throughput sequencing"[DataSet Type] AND ("Homo sapiens"[Organism] OR "Mus musculus"[Organism] OR "rattus norvegicus"[Organism]) AND ("2020/01/01"[PDAT] : "3000"[PDAT])'
$ esearch -db gds -query "$geo_query" | efetch -format uid > gds_results.txt
$ wc -l gds_results.txt
9981 gds_results.txt
$ head -n2 gds_results.txt
200134092
200120931

```

Instead, if you are looking for a way to get the GSE accessions, you can use the built-in `xtract` command to parse the XML returned by `esummary` as follows:

```
$ esearch -db gds -query "$geo_query" | esummary | xtract -pattern DocumentSummary -first Accession > gse_accs.txt > gse_accs.txt
$ wc -l gse_accs.txt
9981 gse_accs.txt
$ head -n2 gse_accs.txt
GSE165829
GSE165824

```

# Answer
> 0 votes
I think I would just try to do this with GEOquery.

https://bioconductor.org/packages/release/bioc/html/GEOquery.html

---
Tags: phylogenetics, geo
---


# Question
Title: What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)?
I bought a "Whole Genome Sequence" kit, which gave me the following VCF files (this service did not include any analysis, just the raw data):

`cnv.vcf` `indel.vcf` `snp.vcf` `sv.vcf`

My understanding is that these files show the differences of my DNA relative to a reference human genome (hg19 in my case)

What tools and resources can I use to find out which known alleles/variants of genes I have?

## An example: Blood type

According to Wikipedia, there are 6 common alleles of the ABO gene, and which alleles I have determine my blood type.

Looking at my VCF results, I have these SNPs and INDELs within that gene:

```
<removed>

```

Given this, how would I figure out my blood type?

I understand that my data may not be 100% accurate, and that many phenotypes are the result of more than one gene, environmental factors, or are not well studied at all, but this should work for some simple things ie blood type, COMT, BHLHE41, CYP1A2, etc

# Answer
> 1 votes
If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.

What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.

There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.

---
Tags: genome, vcf, public-databases, wgs, human-genome
---


title: Fitting branch models with godon

text:
# Question
Title: Fitting branch models with godon
Is it possible to fit branch models (Yang and Nielsen, 1998) in godon? That is, the model which averages over all sites, as opposed to the branch site models?

# Answer
> 2 votes
Branch-models are currently not implemented in godon. In the PAML manual Ziheng Yang writes:

> This model is very parameter-rich and its use is discouraged.

This was one of the reasons I did not implement it.

*Disclaimer: I am the author of godon.*

---
Tags: phylogenetics, codon, godon
---


# Question
Title: Kallisto psudoaligner quant not giving me expected output
I am trying to run this kallisto (latest) quant and all I get are abundance and json files.

```
kallisto quant --plaintext -t 4  -o rf-stranded.nob.1K1_L7.0 -i mm9_ref.idx  --rf-stranded 1K1_L7.0 2K1_L7.0

```

the results files in stderr is:

```
[quant] fragment length distribution will be estimated from the data
[index] k-mer length: 15
[index] number of targets: 23
[index] number of k-mers: 363,704,628
[index] number of equivalence classes: 1,945,936
[quant] running in paired-end mode
[quant] will process pair 1: 1K1_L1.0
                             2K1_L1.0
[quant] finding pseudoalignments for the reads ... done
[quant] processed 22,897,424 reads, 4,476,549 reads pseudoaligned
[quant] estimated average fragment length: 24.2427
[   em] quantifying the abundances ... done
[   em] the Expectation-Maximization algorithm ran for 52 rounds

```

---

I get an abundance file with:

```
target_id       length  eff_length      est_counts      tpm
chr1    197195432       1.89745e+08     253015  21.3249
chr10   129993255       1.31621e+08     210320  25.5544
chr11   121843856       1.38776e+08     263606  30.3773
chr12   121257530       1.23195e+08     143568  18.6371
chr13   120284312       1.22484e+08     132969  17.3613
chr14   125194864       1.20926e+08     112711  14.9058
chr15   103494974       1.07866e+08     109335  16.2102
chr16   98319150        9.70075e+07     109045  17.9767
chr17   95272651        1.06375e+08     169166  25.4322
chr18   90772031        9.23542e+07     121990  21.1241
chr19   61342430        7.23377e+07     157142  34.7407
chr2    181748087       1.81791e+08     244870  21.5414
chr3    159599783       1.45978e+08     209257  22.9247
chr4    155630120       1.61978e+08     226893  22.4014
chr5    152537259       1.65672e+08     194951  18.8186
chr6    149517037       1.43243e+08     246652  27.5373
chr7    152524553       1.82795e+08     197527  17.2812
chr8    131738871       1.50556e+08     203954  21.6643
chr9    124076172       1.30446e+08     222639  27.295
chrX    166650296       1.40278e+08     176259  20.0942
chrY    15902555        5.76797e+07     123.462 0.0342311
mtDNA   16299   11804.8 737285  998819
rDNA    249368  720854  33272.5 738.159

```

run\_info.json:

```
        "n_targets": 23,
        "n_bootstraps": 100,
        "n_processed": 22897424,
        "n_pseudoaligned": 4476549,
        "n_unique": 3641646,
        "p_pseudoaligned": 19.6,
        "p_unique": 15.9,
        "kallisto_version": "0.46.1",
        "index_version": 10,
        "start_time": "",
        "call": "kallisto quant --plaintext -t 4 -b 100 -o rf-stranded.1K1_L1.0 -i mm9_ref.idx --bias --rf-stranded --gtf UCSC_genes.gtf --pseudobam --genomebam 1K1_L1.0.gz 2K1_L1.0.gz"
}

```

---

When I ran it with boot strapping (-b), -gtf and --pseudobam --genomebam it gave me a bam file.

```
44879990 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
8078785 + 0 mapped (18.00% : N/A)
44879990 + 0 paired in sequencing
22439995 + 0 read1
22439995 + 0 read2
8077160 + 0 properly paired (18.00% : N/A)
8077160 + 0 with itself and mate mapped
1625 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

```

I'm not sure what is going wrong, can someone help?

NOTE I use a 15nt kmer, and this about 1/2 of what is suggested (31 I think), is that a problem?

# Answer
As you surmised, kallisto requires a transcriptome fasta file, not a genome fasta file.
> 3 votes

---
Tags: rna-seq, kallisto
---


# Question
Title: Does anyone know how I can convert DNA code into FASTA for this TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC
TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC this is the DNA code Sequence. I have tried converting each into the amino acid and ended up with this LKHWMNEKPCFATPQX but I was told that this isn't right.

# Answer
> 1 votes
This is your nuc seq as a fasta:

```
>sequence_1
TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC

```

To translate it use an online tool (ie https://web.expasy.org/translate/) or do it in biopython:

```
from Bio import SeqIO
for record in SeqIO.parse("/path/seq1.fa", "fasta"):
    print ('nuc', record, sep='\n')
    print ('prot', record.translate(), sep='\n')

```

Which gives:

```
nuc
ID: sequence_1
Name: sequence_1
Description: sequence_1
Number of features: 0
Seq('TTGAAACACTGGATGAATGAAAAGCCCTGCTTTGCAACCCCTCAGC', SingleLetterAlphabet())
prot
ID: <unknown id>
Name: <unknown name>
Description: <unknown description>
Number of features: 0
Seq('LKHWMNEKPCFATPQ', ExtendedIUPACProtein())

```

Note, this defaults to the human codon table. Other codon tables are available with the parameter 'table'.

---
Tags: sequence-alignment, fasta, sequence-analysis, phylogenetics
---


# Question
Title: What are the applications of DNA or RNA pattern matching?
I'm assuming we don't do pattern matching in DNA or RNA for the fun of it. So I'd like to know what are the applications of pattern matching or where does it fit in in larger applications?

I'm a computer scientist and when studying certain topics in CS we see a lot of references to applying algorithms to areas of biology or genetics one of which is applying text pattern matching algorithms to DNA to find sub sequences. Whenever I see those kinds of examples I think "Oh that's pretty cool" but then I don't see how it's applied in the real world for benefit and I can't find anything to suggest where these algorithms are applied in these fields.

# Answer
One useful application of pattern matching is for DNA sequence motifs. As an example of a direct application, there is a class of proteins called transcription factors (TFs), which bind to DNA and regulate the expression of genes. There are many TFs. Each TF likes to bind to its particular conserved sequence of DNA, with some wiggle room for some parts of the sequence.

Algorithms for finding and scoring these patterns (and mismatches in patterns) is important for identifying regions where regulation occurs, which proteins are involved, and for answering biological questions related to gene regulation and basic biological functionality — and malfunction, in case of disease.
> 2 votes

# Answer
These algorithms are used quite a lot, and I will touch only a few aspects:

* **Genetic sequence alignment.** Different (but still closely related) organisms have very similar genetic sequences, which differ only at some positions due to *point mutations*, which can be nucleotide substitutions (another letter in a sequence) or insertion/deletion of a few nucleotides. Many biological questions such as studying new coronavirus strains, predicting genetic diseases, and identifying cancerous cells are based on such sequence comparisons. This requires aligning many sequences, using the appropriate algorithms. This may be not hard for two sequences, but for hundreds or thousands of them it may require serious computational effort.
* **Sequencing reads mapping** Much of the current progress in genetics is due to the advent of the *next generation sequencing*, also called *massively parallel sequencing*, as opposed to Sanger sequencing used in the first human genome project (see here, for example). The basis of this technology is splitting the genetic sequence (a few billion nucleotides in a human) into small fragments of a few hundred nucleotides long and sequencing them simultaneously. One then faces the problem of putting them back together, which can be achieved either via matching them to a known sequence for a similar organism (*sequence mapping*) or by putting the reads together on the basis of their similarity (*sequence assembly*). There exist multiple sophisticated bioinformatic tools for achieving this (BWA, bowtie2, TopHat, Megahit to name just a few, so that you can google them).
* **Blasting** Another application, closely related to my first point, is quickly identifying the sequence as belonging to a certain gene of a certain organism by comparing it to a large database. BLAST is one of the well-known collections of tools for doing this.
* Just to mention (only a bit) more exotic applications: the algorithms associated with parsing context free grammars are used for **RNA structure prediction**.

Durbin's book is a good casual reading (for somebody with a computer science background) to get an overview without knowing too much biology.
> 2 votes

# Answer
The applications can be categorized based on whether you study the sequence for analysis, or for DNA synthesis and assembly.

In addition to TF patterns mentioned in this answer, other patterns for DNA binding proteins, enchancers, promoters, and the Kozak sequence (`GCCGCC[A|G]CCATGG..`) may also be useful for discovery of ORFs (genes) / regulation.

Certain methylase enzymes recognize a DNA pattern. For example EcoKI methylates adenine in `AACNNNNNNGTGC` (2nd A) and its reverse complement, `GCACNNNNNNGTT` (3rd A). Another such enzyme is EcoBI. This is useful for finding methylation sites, which may interfere with DNA restriction or expression.

Some restriction enzymes also recognize a DNA pattern: HinFI recognizes

```
5'GANTC
3'CTNAG

```

Other such enzymes are EcoP15I, EcoRII. Finding these sites are very important for designing a DNA cloning or manipulation experiment.

For DNA synthesis, it is useful to avoid homopolymers (`N+`) and repeated k-mers and hairpins (a sequence segment which has a reverse complement "nearby") as these can cause problems during synthesis.

Direct, tandem or inverted repeats can cause problems during cloning.

Finally, a perhaps trivial example is matching a polypeptide sequence, using DNA, as the genetic code is redundant. An example would be searching for target peptides in a protein.
> 1 votes

# Answer
To extend on the existing answers I would say that pattern matching algorithms are at the core of bioinformatics, certainly any NGS realted bioinformatics. The Smith-Waterman algorithm and its variants are probably the most used https://en.wikipedia.org/wiki/Smith%E2%80%93Waterman\_algorithm .

Reference genome assembly (the human genome project for example) is based on kmer based sequence matching. Variant calling is based on mapping sequence to this reference, so thats pretty much all clinical bioinformatics. And a million other tools - one random example is dectecting horizontal gene transfer of antibiotic resistence genes.
> 1 votes

---
Tags: rna-seq, phylogenetics
---


# Question
Title: How to convert a SDF molecular structure file to a pandas dataframe with Python?
What is the best python tool to convert a SDF file to a structured data frame (pandas.DataFrame). Similar to Molconvert from ChemAxon? SDF can have different datafields for each compound which I would like to have in a dataframe. The 2D or 3D structures can be disregarded.

# Answer
> 2 votes
It is easily done with RDKit:

```
from rdkit.Chem import PandasTools

fn = '/mypath/dataset.sdf'

df = PandasTools.LoadSDF(fn, embedProps=True, molColName=None, smilesName='smiles')

```

---
Tags: python, rdkit
---


title: blast p-error with making directory

text:
# Question
Title: blast p-error with making directory
I tried running this command `python blastp.py CmydPro tlr TLR.txt` to extract the query sequence for the target species I'm looking at. However, whenever I try running it always comes up as a blast p error and says there's no directory. I downloaded blastp.py as a file made from my lab's GitHub but I have no idea how to troubleshoot the error since it's not a required directory. I made sure all my directories were correct and that all my files were intact. Any help with this would be greatly appreciated.

This is the error I keep getting:

```
Traceback (most recent call last):
File "blastp.py", line 53, in
blastp_scout = subprocess.Popen(['blastp', '-db', organism_db, '-query', 'temp.fa', '-out', 'blast', '-num_descriptions', '5', '-num_alignments', '0'])
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 854, in __init__
self._execute_child(args, executable, preexec_fn, close_fds,
File "/home/dylan_w190/miniconda3/lib/python3.8/subprocess.py", line 1702, in _execute_child
raise child_exception_type(errno_num, err_msg, err_filename)
FileNotFoundError: [Errno 2] No such file or directory: 'blastp'

```

These are all the files I have stored in my directory:

```
'$RECYCLE.BIN' 'Genome Files' TLR.txt
'ASM 28112' 'Genome Files (1)' Tham
'ASM 28964' 'Genome Files (2)' Vim
'ASM 3597' 'Genome Files (3)' Xenopus
'ASM 45574' 'Genome Files (4)' blast-2.11.0+
AnoCar Git blastp.py
Chrysemys GopAga get-pip.py
Cliv GopAga2.0-Genome human
CmydPro Halia ls
CmydPro_results.txt LICENCE mouse
Config.Msi Lib ncbi-blast-2.6.0+
CroPor Miniconda3-latest-Linux-x86_64.sh ncbi-blast-2.6.0+-win64.exe
'FASTAs (2)' OpenSSH-Win64 pip
Format Ornitho pvi
GCF_015237465.1_rCheMyd1.pri_genomic.fna PelSin pyzo
GFFs Python setup.py
'GFFs (2)' PythonX temp.fa
GRCh38 README.txt tlr_CmydPro.fa
GRCm38 Rnor tlr_ref.fasta
Gallus 'Rnor - Copy' vs_buildtools__576366542.1609221376.exe
GavGan Scripts vs_community__576366542.1609221376.exe
Gekko 'System Volume Information'

```

There is a windows version. It seems that I already downloaded blast but I need to do this now:

`Command 'blastp' not found, but can be installed with: sudo apt install ncbi-blast+`

Do you have any idea how to might download this?

I already downloaded blast as a tar file and the 2.11.0+ version. I'm just not sure how to download just ncbi-blast+ on its own.

# Answer
You don't have blastp in your path. Blast is either not installed on your system or it's just not been added to path right.
> 1 votes

---
Tags: file-formats
---


# Question
Title: Is there a database of protein sequences/structure along with their melting temperature?
Such databases have been constructed for example in here or here. But I can't find them anywhere online. Can someone point me to where I can find one available online?

Otherwise I'll have to write to these people and I don't know when/whether I'll get a reply.

# Answer
> 1 votes
ProTherm2 contains thermodynamic properties of proteins and mutants. https://www.iitm.ac.in/bioinfo/ProTherm2/ .

NAR Database issue advertises it as containing melting temperatures.

For the reference - found it with https://biokeanos.com/search?q=protein+%2Bmelting

---
Tags: database, public-databases
---


# Question
Title: Plotting gene models - alternatives to ggbio
I am curious how people make genome track plots such as what can be done with the Bioconductor package `ggbio`.

For example this plot was made w/ ggbio:

I think ggbio is great for making fast plots, but it is kind of a pain to make publication quality figures because it is hard to manipulate the plots after they are constructed. I have noticed that ggbio constructed plots do not return a ggplot or gtable object. I can access the grobs from a ggbio tracks object, but they won't be aligned as they would be in the final plot. I like to make shared legends and include them in the plot, and without having a gtable or ggplot object it is hard to add the shared legends in, and I end up having to save them as a separate file which is tedious. FYI I mainly use ggbio for gene tracks + x-axis alignment.

Normally I would use ggplot2 to customize the plot extensively, but it doesn't work very well after plotting with ggbio.

How can I use ggbio for these types of plots? Is there a better tool instead?

EDIT: I would get the legends for each individual plot and combine them before adding back to the main plot using something like this:

```
leg1.grob = cowplot::get_legend(plot1)
leg2.grob = cowplot::get_legend(plot2)
sharedleg.grob = gridExtra::gtable_rbind(leg1,leg2)

mainplot.grob = ggplotGrob(main_plot)
mainplot.grob = gtable_add_cols(mainplot.grob,
                                widths = sum(sharedleg.grob$widths))
mainplot.grob = gtable_add_grobs(mainplot.grob,
                                 grobs = sharedleg.grob,
                                 t=t,b=b,l=l,r=r)
cowplot::ggdraw(mainplot.grob)

```

I've used gviz before, but long ago. Do you know if you can extract the grid plot grobs of the gene track plots with Gviz?

# Answer
> 1 votes
It seems like the ggbio object already contains the corresponding ggplot object. You should be able to extract the ggplot from the `ggbio_obj@ggplot` slot and extract your gtables.

Small example:

```
library(ggbio)
library(ggplot2)
library(patchwork)

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- as(c("chr17:41239394-41319151:+"), "GRanges")
p_cov <- autoplot(fl.bam, which = wh)
#> reading in as Bamfile
#> Parsing raw coverage...
#> Read GAlignments from BamFile...
#> extracting information...

df <-  data.frame(x= sample(1:100,size = 100,replace = TRUE))
p_hist <- ggplot(df) +
    geom_histogram(aes(x))

p_hist / p_cov
#> Error: Can't add `e2` to a ggplot object.
#> Backtrace:
#>     █
#>  1. └─patchwork:::`/.ggplot`(p_hist, p_cov)
#>  2.   └─GGally:::`+.gg`(e1, e2)
#>  3.     └─e1 %+% e2
#>  4.       └─ggplot2:::add_ggplot(e1, e2, e2name)
#>  5.         ├─ggplot2::ggplot_add(object, p, objectname)
#>  6.         └─ggplot2:::ggplot_add.default(object, p, objectname)

p_hist / p_cov@ggplot
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

```

Created on 2021-02-05 by the reprex package (v1.0.0)

Session info
```
sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Red Hat Enterprise Linux
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Berlin
#>  date     2021-02-05
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version  date       lib source
#>  acepack                1.4.1    2016-10-29 [1] CRAN (R 3.6.0)
#>  AnnotationDbi          1.48.0   2019-10-29 [1] Bioconductor
#>  AnnotationFilter       1.10.0   2019-10-29 [1] Bioconductor
#>  askpass                1.1      2019-01-13 [1] CRAN (R 3.6.0)
#>  assertthat             0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
#>  backports              1.1.8    2020-06-17 [1] CRAN (R 3.6.0)
#>  base64enc              0.1-3    2015-07-28 [1] CRAN (R 3.6.0)
#>  Biobase                2.46.0   2019-10-29 [1] Bioconductor
#>  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor
#>  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor
#>  BiocManager            1.30.10  2019-11-16 [1] CRAN (R 3.6.0)
#>  BiocParallel           1.20.1   2019-12-21 [1] Bioconductor
#>  biomaRt                2.42.1   2020-03-26 [1] Bioconductor
#>  Biostrings             2.54.0   2019-10-29 [1] Bioconductor
#>  biovizBase             1.34.1   2019-12-04 [1] Bioconductor
#>  bit                    1.1-15.2 2020-02-10 [1] CRAN (R 3.6.0)
#>  bit64                  0.9-7    2017-05-08 [1] CRAN (R 3.6.0)
#>  bitops                 1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
#>  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.0)
#>  BSgenome               1.54.0   2019-10-29 [1] Bioconductor
#>  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.0)
#>  cli                    2.0.2    2020-02-28 [1] CRAN (R 3.6.0)
#>  cluster                2.0.8    2019-04-05 [2] CRAN (R 3.6.0)
#>  colorspace             1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
#>  crayon                 1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
#>  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.0)
#>  data.table             1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
#>  DBI                    1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
#>  dbplyr                 1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
#>  DelayedArray           0.12.2   2020-01-06 [1] Bioconductor
#>  dichromat              2.0-0    2013-01-24 [1] CRAN (R 3.6.0)
#>  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.0)
#>  dplyr                  0.8.5    2020-03-07 [1] CRAN (R 3.6.0)
#>  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.0)
#>  ensembldb              2.10.2   2019-11-20 [1] Bioconductor
#>  evaluate               0.14     2019-05-28 [1] CRAN (R 3.6.0)
#>  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.0)
#>  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.0)
#>  foreign                0.8-71   2018-07-20 [2] CRAN (R 3.6.0)
#>  Formula                1.2-3    2018-05-03 [1] CRAN (R 3.6.0)
#>  fs                     1.4.0    2020-03-31 [1] CRAN (R 3.6.0)
#>  GenomeInfoDb           1.22.1   2020-03-27 [1] Bioconductor
#>  GenomeInfoDbData       1.2.2    2020-04-01 [1] Bioconductor
#>  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor
#>  GenomicFeatures        1.38.2   2020-02-15 [1] Bioconductor
#>  GenomicRanges          1.38.0   2019-10-29 [1] Bioconductor
#>  GGally                 1.5.0    2020-03-25 [1] CRAN (R 3.6.0)
#>  ggbio                * 1.34.0   2019-10-29 [1] Bioconductor
#>  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.0)
#>  glue                   1.4.1    2020-05-13 [1] CRAN (R 3.6.0)
#>  graph                  1.64.0   2019-10-29 [1] Bioconductor
#>  gridExtra              2.3      2017-09-09 [1] CRAN (R 3.6.0)
#>  gtable                 0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
#>  highr                  0.8      2019-03-20 [1] CRAN (R 3.6.0)
#>  Hmisc                  4.4-0    2020-03-23 [1] CRAN (R 3.6.0)
#>  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.0)
#>  htmlTable              1.13.3   2019-12-04 [1] CRAN (R 3.6.0)
#>  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.0)
#>  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.0)
#>  httr                   1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
#>  IRanges                2.20.2   2020-01-13 [1] Bioconductor
#>  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.0)
#>  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.0)
#>  labeling               0.3      2014-08-23 [1] CRAN (R 3.6.0)
#>  lattice                0.20-38  2018-11-04 [2] CRAN (R 3.6.0)
#>  latticeExtra           0.6-29   2019-12-19 [1] CRAN (R 3.6.0)
#>  lazyeval               0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
#>  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.0)
#>  magrittr               1.5      2014-11-22 [1] CRAN (R 3.6.0)
#>  Matrix                 1.2-17   2019-03-22 [2] CRAN (R 3.6.0)
#>  matrixStats            0.56.0   2020-03-13 [1] CRAN (R 3.6.0)
#>  memoise                1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
#>  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.0)
#>  munsell                0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
#>  nnet                   7.3-12   2016-02-02 [2] CRAN (R 3.6.0)
#>  openssl                1.4.1    2019-07-18 [1] CRAN (R 3.6.0)
#>  OrganismDbi            1.28.0   2019-10-29 [1] Bioconductor
#>  patchwork            * 1.0.0    2019-12-01 [1] CRAN (R 3.6.0)
#>  pillar                 1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
#>  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
#>  plyr                   1.8.6    2020-03-03 [1] CRAN (R 3.6.0)
#>  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.0)
#>  progress               1.2.2    2019-05-16 [1] CRAN (R 3.6.0)
#>  ProtGenerics           1.18.0   2019-10-29 [1] Bioconductor
#>  purrr                  0.3.4    2020-04-17 [1] CRAN (R 3.6.0)
#>  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
#>  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.0)
#>  RBGL                   1.62.1   2019-10-30 [1] Bioconductor
#>  RColorBrewer           1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp                   1.0.4.6  2020-04-09 [1] CRAN (R 3.6.0)
#>  RCurl                  1.98-1.1 2020-01-19 [1] CRAN (R 3.6.0)
#>  reprex                 1.0.0    2021-01-27 [1] CRAN (R 3.6.0)
#>  reshape                0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
#>  reshape2               1.4.4    2020-04-09 [1] CRAN (R 3.6.0)
#>  rlang                  0.4.6    2020-05-02 [1] CRAN (R 3.6.0)
#>  rmarkdown              2.1      2020-01-20 [1] CRAN (R 3.6.0)
#>  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.0)
#>  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor
#>  RSQLite                2.2.0    2020-01-07 [1] CRAN (R 3.6.0)
#>  rstudioapi             0.11     2020-02-07 [1] CRAN (R 3.6.0)
#>  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor
#>  S4Vectors              0.24.3   2020-01-18 [1] Bioconductor
#>  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.0)
#>  sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi                1.4.6    2020-02-17 [1] CRAN (R 3.6.0)
#>  stringr                1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
#>  styler                 1.3.2    2020-02-23 [1] CRAN (R 3.6.0)
#>  SummarizedExperiment   1.16.1   2019-12-19 [1] Bioconductor
#>  survival               3.1-11   2020-03-07 [1] CRAN (R 3.6.0)
#>  tibble                 3.0.0    2020-03-30 [1] CRAN (R 3.6.0)
#>  tidyselect             1.0.0    2020-01-27 [1] CRAN (R 3.6.0)
#>  VariantAnnotation      1.32.0   2019-10-29 [1] Bioconductor
#>  vctrs                  0.2.4    2020-03-10 [1] CRAN (R 3.6.0)
#>  withr                  2.4.1    2021-01-26 [1] CRAN (R 3.6.0)
#>  xfun                   0.15     2020-06-21 [1] CRAN (R 3.6.0)
#>  XML                    3.99-0.3 2020-01-20 [1] CRAN (R 3.6.0)
#>  xml2                   1.3.0    2020-04-01 [1] CRAN (R 3.6.0)
#>  XVector                0.26.0   2019-10-29 [1] Bioconductor
#>  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.0)
#>  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor

```

---
Tags: r, visualization, ggplot2
---


# Question
Title: A good tool for gene locus visualization
Does anyone know of a good tool for visualizing gene loci? Basically, I'd like to show how a group of genes look in different strains of a bacteria. I have annotation of course, I just need a tool to plug it into and get an image of the locus to compare.

I'm talking something basic like this:

# Answer
> 3 votes
I ended up using Easyfig, which is made for the task and requires little tinkering around, once fed the genbank files. (Granted, being so simple, it's hard to make the image more advanced, but it serves my purposes for now).

Here's an example of my image.

# Answer
> 5 votes
I think in R, Gviz is the best choice. I have used it for human sequences but it will probably also be suitable for Prokaryotes.

Another option is MS powerpoint, if you only need the arrows...

# Answer
> 3 votes
GFF2PS was developed for exactly this sort of thing. It takes a GFF file as input, and will produce a pretty PS output showing the annotations. It was used to create the poster for Celera's human genome draft \[1\].

Here's a very simple example of what it can do:

And a more complex one from the mosquito genome paper \[2\] (click here for a full size, 1.5M gzipped PostScript file):

If you are happy with the defaults, it is very easy to use, but otherwise it is a bit of a pain to configure. It is absurdly powerful though.

---

1. Venter et al. Science 291(5507):1304-1351 (2001)
2. Holt et al. Science 298(5591):129-149 (2002)

# Answer
> 3 votes
gggenes is a ggplot2 based solution and looks very appealing. Plus, it has a nice introduction with examples.

```
library(ggplot2)
library(gggenes)
ggplot(example_genes, aes(xmin = start, xmax = end, y = molecule, fill = gene)) +
  geom_gene_arrow() +
  facet_wrap(~ molecule, scales = "free", ncol = 1) +
  scale_fill_brewer(palette = "Set3") +
  theme_genes()

```

---
Tags: gene, visualization
---


# Question
Title: How to remove sequences from a fasta file using a sequence ID list which contains a space within the id?
I have a fasta file that contains sequence reads and sequence id file that needed to be removed from the fasta file. I have done this earlier, but since id contains a space my piece of code is not working and I don't know how to change it since I'm still learning. fasta file looks like this.

```
>m64041_200717_231916/74/12723_24868 id=30
CCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCA
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/105/20691_65844 id=488
AGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCAAGACTTTAT
>m64041_200717_231916/108/17414_67048 id=4956
TGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

```

IDList file is like this.

```
m64041_200717_231916/74/12723_24868 id=30
m64041_200717_231916/108/17414_67048 id=4956
m64041_200717_231916/105/20691_65844 id=488

```

I need to remove these sequences from that fasta file. I have a code I used earlier. But it is not working since id contains a space (I think). Previous code is like this.

```
# use as follows
# ./filter_seq.bash <fastafile> <list_of_ids_you_want_to_remove>
# ./filter_seq.bash all.fasta bacterial_IDlist

grep ">" $1 | sed 's/>//g' | sort | uniq > all_del
sort $2 | uniq > query_del
comm -23 all_del query_del | sort | uniq > filtered_del
perl -ne 'if(/^>(\S+)/){$c=$i{$1}}$c?print:chomp;$i{$_}=1 if @ARGV' filtered_del $1 > filtered.fasta
rm *_del

```

I would be grateful if someone can help me.

# Answer
> 1 votes
Install (pip) biopython if u don't have it.

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

```

# Answer
> 0 votes
I tend to do this using the FastaToTbl and TblToFasta scripts I have posted before. `FastaToTbl` will convert fastq sequences to single-line, tab-separated recors which makes them far easier to parse. Using these two scripts, you can do:

```
$ FastaToTbl file.fa | grep -vwf idList | TblToFasta
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCC
TGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCAC
ACTTG

```

The `-v` option will tell `grep` to return non-matching lines, and the `-f` option tells it to read the patterns to search for from a file. In this case, we give it the file with the IDs to skip (`idList`). The `-w` ensures we only find "whole word" matches, so that something like `m64041_200717_231916/77/1941_50622 id=31151111` will not count as a match for `m64041_200717_231916/77/1941_50622 id=31151111`.

---

Another approach would be to just use `awk` directly:

```
$ awk '(NR==FNR){s[">"$0]++} (/^>/){ s[$0] ? a=0 : a=1}a' idList file.fa
>m64041_200717_231916/77/1941_50622 id=3115
TGAGCCGCACCTCCTCAATCTGCAGCAGTTGAGGCCACTACCCTTCTGCTCAATGGTTCCTGCAGACTTTATCATAGTCACTCACACTTGTCCATGAG
>m64041_200717_231916/162/0_6615 id=857
CAGACTTTATCATAGTCACTCACACTTGTCCATGAGGCAGACTTTATCATAGTCACTCACACTTG

```

Here is the same `awk` script slightly expanded to make it easier to understand:

```
 awk '{
        if (NR==FNR){idsToSkip[">"$0]++}
        else if (/^>/){
            if($0 in idsToSkip){
                wantThisId = 0
            }
            else{
                wantThisId = 1
            }
        }
        if(wantThisId == 1){
            print
        }
      }' idList file.fa

```

---
Tags: bash, linux, filtering
---


title: creating a tab delimited file

text:
# Question
Title: creating a tab delimited file
I am working on a project using a fasta file. I am writing my command in nano within command-line and executing using python, also within my command-line.

I would like my command to provide me with a tab delimited file with three columns: first column should contain my sequence name, second column should provide me with my sequence length, and the third column should show the sequence itself.

I have written the following command so far within nano:

```
from Bio import SeqIO
import sys
for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):

  list = hello_fasta.split("\t")

  print hello_fasta.description
  print (len(hello_fasta.seq))

```

For example, I would like my command to provide me with the desired output and with the following order: Gene name ; Gene length ; Gene seq

H0192X 26 FORUWOHRPPTRWFAWWEAKJNFWEJ

# Answer
> 1 votes
You can use a list and `insert()` to add an element in a specific order, then expand the list with `*`. Or you can use `join()`.

```
from Bio import SeqIO
import sys

for hello_fasta in SeqIO.parse(sys.argv[1], "fasta"):
  sequences = []
  sequences.insert(0, hello_fasta.description)
  sequences.insert(1, len(hello_fasta.seq))
  sequences.insert(2, hello_fasta.seq)
  # option 1
  print(*sequences, sep='\t')
  # option 2
  print('\t'.join(map(str, sequences)))

```

# Answer
> 0 votes
Here's a solution using `pandas` if you want to save the tsv:

```
from Bio import SeqIO
import pandas as pd
from io import StringIO

example = """
>seq0
FQTWEEFSRAAEKLYLADPMKVRVVLKYRHVDGNLCIKVTDDLVCLVYRTDQAQDVKKIEKF
>seq1
KYRTWEEFTRAAEKLYQADPMKVRVVLKYRHCDGNLCIKVTDDVVCLLYRTDQAQDVKKIEKFHSQLMRLMELKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
>seq2
EEYQTWEEFARAAEKLYLTDPMKVRVVLKYRHCDGNLCMKVTDDAVCLQYKTDQAQDVKKVEKLHGK
>seq3
MYQVWEEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVCLQYKTDQAQDVK
>seq4
EEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVVSYEMRLFGVQKDNFALEHSLL
>seq5
SWEEFAKAAEVLYLEDPMKCRMCTKYRHVDHKLVVKLTDNHTVLKYVTDMAQDVKKIEKLTTLLMR
>seq6
FTNWEEFAKAAERLHSANPEKCRFVTKYNHTKGELVLKLTDDVVCLQYSTNQLQDVKKLEKLSSTLLRSI
>seq7
SWEEFVERSVQLFRGDPNATRYVMKYRHCEGKLVLKVTDDRECLKFKTDQAQDAKKMEKLNNIFF
>seq8
SWDEFVDRSVQLFRADPESTRYVMKYRHCDGKLVLKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
>seq9
KNWEDFEIAAENMYMANPQNCRYTMKYVHSKGHILLKMSDNVKCVQYRAENMPDLKK
>seq10
FDSWDEFVSKSVELFRNHPDTTRYVVKYRHCEGKLVLKVTDNHECLKFKTDQAQDAKKMEK
"""

# This example just happens to be a string, just load your
# fasta file using the method you're already using
example_records = SeqIO.parse( StringIO(example), 'fasta')

# Dictionary to hold the data you eventually want in the tsv
data = {"Gene name" : list(),
        "Gene length" : list(),
        "Gene seq" : list()}

# Append the necessary into the data dictionary
for record in example_records:
    data['Gene name'].append(record.description)
    data['Gene length'].append(len(record.seq))
    data['Gene seq'].append(str(record.seq))

# Convert your data into a pandas DataFrame and save as a tsv
gene_df = pd.DataFrame(data)
gene_df.to_csv("gene_info.tsv", sep = '\t', index = False)

```

This results in a tsv that looks like this:

```
$ head gene_info.tsv
Gene name       Gene length     Gene seq
seq0    62      FQTWEEFSRAAEKLYLADPMKVRVVLKYRHVDGNLCIKVTDDLVCLVYRTDQAQDVKKIEKF
seq1    106     KYRTWEEFTRAAEKLYQADPMKVRVVLKYRHCDGNLCIKVTDDVVCLLYRTDQAQDVKKIEKFHSQLMRLMELKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM
seq2    67      EEYQTWEEFARAAEKLYLTDPMKVRVVLKYRHCDGNLCMKVTDDAVCLQYKTDQAQDVKKVEKLHGK
seq3    58      MYQVWEEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVCLQYKTDQAQDVK
seq4    62      EEFSRAVEKLYLTDPMKVRVVLKYRHCDGNLCIKVTDNSVVSYEMRLFGVQKDNFALEHSLL
seq5    66      SWEEFAKAAEVLYLEDPMKCRMCTKYRHVDHKLVVKLTDNHTVLKYVTDMAQDVKKIEKLTTLLMR
seq6    70      FTNWEEFAKAAERLHSANPEKCRFVTKYNHTKGELVLKLTDDVVCLQYSTNQLQDVKKLEKLSSTLLRSI
seq7    65      SWEEFVERSVQLFRGDPNATRYVMKYRHCEGKLVLKVTDDRECLKFKTDQAQDAKKMEKLNNIFF
seq8    68      SWDEFVDRSVQLFRADPESTRYVMKYRHCDGKLVLKVTDNKECLKFKTDQAQEAKKMEKLNNIFFTLM

```

Hopefully this helps!

---
Tags: python, assembly, file-formats
---


title: Diploid Consensus Sequences

text:
# Question
Title: Diploid Consensus Sequences
My goal is to generate a diploid consensus genome. By "diploid consensus", I mean that I want to merge all supporting reads for a genomic region into a single haplome such that the resulting sam file contains TWO sets of reads, one from each haplome.

I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species.

I am worried that, by simply using Picard's MarkDuplicates with the "REMOVE\_DUPLICATES" option set to "TRUE", I will retain only the homolog with the highest quality and output a haploid version of my genome. Here is my current workflow:

```

## Map Illumina reads to the reference genome. Include read IDs
bwa mem -M -R '@RG\tID:BioSample.LibraryName\tSM:BioSample\tLB:LibraryName\tPU:SingleFlowcell.SingleLane\tPL:ILLUMINA\t' reference.fasta SRR[number]_1.fastq SRR[number]_2.fastq > SRR[number].sam

## Sort Illumina reads by coordinate
java -jar $EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR

## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar $EBROOTPICARD/picard.jar MarkDuplicates REMOVE_DUPLICATES=true REMOVE_SEQUENCING_DUPLICATES=true I= 'sorted_SRR[number].sam' O= nodup_SRR[number].sam METRICS_FILE= metrics_SRR[number].txt TMP_DIR=$TMPDIR

## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam > nodup_SRR[number].bam

## Re-align reads around indels to reduce slippage.
gatk LeftAlignIndels -R reference.fasta -I nodup_SRR[number].bam -O aligned_SRR[number].bam

```

Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?

# Answer
> 1 votes
How to separate reads, one for each reference haplotype.

1. Map all the reads against the reference genome
2. Use GATK to call SNPs
3. Use GATK tool \[FastaAlternateReferenceMaker\]\[1\] to build the alternative reference
4. You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
5. Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called \[BayesASE\]\[2\]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved. \[1\]: https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker \[2\]: https://github.com/McIntyre-Lab/BayesASE

---
Tags: sam, markduplicate
---


# Question
Title: Error in if (is.na(sra_acc$run[i])) { : argument is of length zero
I am trying to write a code for importing sra files as fastq files in Rstudio (for obtaining and aligning RNA-seq reads from public repositories). I keep getting "Error in if (is.na(sra\_acc$run\[i\])) { : argument is of length zero" right after I run this code. Any suggestions on what might be going wrong? My code is as follows:

```
library(SRAdb)
sqlfile <-'SRAmetadb.sqlite'
if(!file.exists('SRAmetadb.sqlite')) sqlfile <<- getSRAdbFile()
sra_con <- dbConnect(SQLite(),sqlfile)
sraInf <- getSRAinfo("SRP045534",sra_con, sraType="sra")

Error in if (is.na(sra_acc$run[i])) { : argument is of length zero

```

# Answer
> 2 votes
After checking with `file.exists()`, I would suggest you try the same command with another SRA run. It may happen that some runs are not accessible for reasons I do not know. I always have problems in accessing sra formatted files and use a workaround to download fastq.

Here is my not official workaround:

```
runlist<-sraConvert("SRP045534",sra_con=sra_con)
runlist[1:2,]

study submission    sample experiment        run
1 SRP045534  SRA179382 SRS685546  SRX681992 SRR1552451
2 SRP045534  SRA179382 SRS685544  SRX681990 SRR1552449

for (dloadrun in runlist$run)
{
    getSRAfile(dloadrun,sra_con=sra_con,fileType="fastq")
}

```

This should iteratively download all the runs associated to the project you are studying. It is possible that you can do the same with sra type file, but I have never been able to use them with sraDB, and never investigated the issue.

---
Tags: r, sratoolkit
---


# Question
Title: Is there a difference between any scRNA-seq visualization and pseudotime?
Is there a good definition of pseudotime? Some tools are clearly labeled as pseudotime and produce values along a trajectory, but there are more complex approaches that involve branching and disconnected populations. There are 70 of them in the dynverse, so there are many options.

UMAP will also order cells, but that is generally not considered pseudotime. You can even overlay Monocle or PAGA results onto a UMAP so the cells are still ordered the same way. SPRING is specifically designed to capture continuous trajectories, but is also not considered pseudotime. Why not?

# Answer
> 1 votes
Pseudotime is the order of cells along a trajectory which in the simplest case is a minimum-spanning tree based on a (low-dimensional) manifold such as PCA or UMAP. Without a defined trajectory you have no pseudotime, and therefore UMAP is not pseudotime, but in fact can be used to define the trajectory.

---
Tags: scrnaseq, visualization, umap, monocle
---


# Question
Title: Retrieving all genes for a gene ontology term
I'm looking for an easy way to retrieve all the genes in a list that are associated with a certain GO term, preferably using R/Bioconductor packages. I'm not interested in under/overrepresentation or enrichment.

For instance, I want a list of all genes known to be located in 'presynaptic endosome' (GO:009883).

I tried the method referred to in an older post (Retrieve all genes for a specific gene ontology and https://www.biostars.org/p/52101/); following is my code:

```
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") #retrieve human.ensembl data via biomaRt

go_1 = c("GO:0098830", "GO:0098954", "GO:0099007", "GO:0099067", "GO:0099037", "GO:0098955", "GO:0099592", "GO:0099532") #GO IDs of all terms associated with "presynaptic endosome" (inclusive of child terms)

pre.gene.data <- getBM(attributes=c('hgnc_symbol', 'ensembl_gene_id', 'go_id', 'go_linkage_type'),
                    filters = 'go', values = list(go_1), mart = ensembl)

```

but my code doesn't filter the genes based on the GO terms. It gives me the following output:

```
GO:0006886 AP3D1 ENSG00000065000 intracellular protein transport IEA
GO:0016192 AP3D1 ENSG00000065000 vesicle-mediated transport IEA
GO:0030117 AP3D1 ENSG00000065000 membrane coat IEA

... and 49 other entries for the gene "AP3D1"

```

I'm unable to understand what might be the problem. I replaced the filter "go\_id" to "go" because of the updated filters used:listFilters(ensembl).

Please help!

# Answer
It is filtering by your GO terms, it is only giving you genes associated with those GO terms. However the GO terms in the output are not linked to the GO terms in your input, so it is giving you all the GO terms linked to the genes which are associated with the GO terms you filtered by.

Essentially BioMart is not very clever. It has a table of all the GO terms linked to the genes, it filters the genes to only include those linked to the GO terms you've given it. Then it takes that list of genes and queries that same table completely separately, getting the GO terms for each gene.
> 1 votes

---
Tags: bioconductor, gene, biomart, gene-ontology
---


title: Supplementary aligments in VAF

text:
# Question
Title: Supplementary aligments in VAF
*This question has also been asked on Biostars*

I have a doubt, are supplementary alignment usually considered when the variant allele frequency is calculated? Thanks a lot.

In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I would like to understand if the supplementary reads aligned in that position must be taken into consideration.

I've used pysamstats to calculate and it consider them, but my doubt is more if it right to consider them, if it is normally done or not. We are doing it on human and we expect only a change of a single nucleotide variant in the region if the editing works.

I'm interested in what is done generally, if there is a "rule" to consider or not the supplementary alignment when the AF is calculated. In my case I have some Illumina short data edited with Crisp-Cas9 and I want to calculate for a specific position all the frequency of the 4 nucleotides. I was trying to understand if the supplementary reads aligned in that position must be taken into consideration.

GATK HaplotypeCaller and Mutect2 take them into consideration? I'm not looking for a specific variant caller, but if I want to take one variant caller as "reference", I would probably look at them.

Thanks a lot.

# Answer
> 1 votes
It depends on the tool and the use case but generally they should be.

It's not as important with short reads as only a small portion of reads have supplementary mapping (unless you're looking at breast cancer post chromothripsis, for example). But imagine the case where you're working with long reads (as I do). I have samples with a mean read length of 60kb. The human reference is structurally different to any given sample meaning that the chance of any given long read having supplementary mapping is high. If supplementary mapping didn't count towards allele counts then many would be artificially low. This is my use case, others are different.

I would count them unless you have some reason to suspect that they are are artefact of the editing and or the mapping. Sorry to never give a straight answer but that's bioinformatics! From what you described I would guess 99% are real but there's normally some noise and if that's of concern your code should include some QC. I.e, consider if multiple reads support the sup mapping and their quality.

---
Tags: sequence-alignment, snp, ngs, variant-calling
---


# Question
Title: Are there open-sourced graph-based variant callers?
I would like to compare GATK with graph-based variant callers. I have seen the Fast and accurate genomic analyses using genome graphs paper from SevenBridges and the paragraph tool by Illumina, though it does not look like folks use their platform or it is not open-sourced (in the case of SevenBridges's GRAF).

The idea is to see if graph-based approaches detect better both SNPs and structural variants compared with existing methods. Additionally, I believe that if there exists some novel graph-based aligner/assembler one could use the metadata downstream to identify variants.

# Answer
Yes, there are. There are some suggestions in the comments (VG, WhatsHap, GraphAligner, Minigraph).

However, to be clear, the current default variant calling algorithm used in the GATK (the haplotype caller) is a graph-based algorithm. It constructs a De Bruijn graph from mapped reads for local assembly around variants. Other tools that also follow this strategy include Platypus and freebayes.

You might also consider multi-color genome assemblers. For example, Cortex and Corticall use colored De Bruijn graphs to aggregate reads from multiple genomes and call variants.
> 4 votes

---
Tags: variant-calling, reference-genome, graphs
---


# Question
Title: Books on application of Machine Learning in Bioinformatics
I'm looking for books on applications of machine learning and statistical data mining in bioinformatics with example codes provided in R and/or python.

I came across the following two books: the first being a quite good one, provides the code in MATLAB, which I'm not very familiar with. The second one too is a good book with examples in R.

1. Feature Selection and Ensemble Methods for Bioinformatics by Oleg Okun
2. Data Science and Predictive Analytics: Biomedical and Health Applications using R by Ivo Dinov

Can people suggest other books based on their first hand experiences?

Thanks

# Answer
> -2 votes
There are plenty of books when we talk about bioinformatics but only a few of them are available for application of machine learning in bioinformatics. From these selected books I would suggest that you go ahead with Statistical Modelling and Machine Learning Principles for Bioinformatics Techniques, Tools, and Applications. Its a recent book published in 2020 and covers most of the important topics in bioinformatics along with their ML applications. Apart from that you can also refer Machine Learning in Bioinformatics it’s quite old but still it covers most of the concepts with in-depth explanation. Along with that you can refer to this amazing book Computational Genomics with R available free online. This will cover both theory and code for different ML algorithms (both supervised and unsupervised) with their application in bioinformatics. It also covers important topics like RNA Seq Analysis, ChIP Seq Analysis in depth with code explanation of different real world examples which gives you head start in this beautiful subject.

# Answer
> 0 votes
One could read *Bioinformatics: The machine learning approach* as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...

...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.

I would like however to stress the importance of the basic background - one regularly sees announcements about *breakthroughs* achieved by reinventing techniques that were already described in the book that I mentioned.

---
Tags: machine-learning
---


# Question
Title: How to perform bowtie2 analysis with slurm?
I am trying to run my alignment script that works locally, using `SBATCH`.

Official manual for `bowtie2` says I can use `samtools` for that. However when I run the following script I have an error

```
line 22: bowtie2: command not found

```

Below is the script:

```
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
module load samtools/1.9

for FILE in *R1_001.fastq.gz;

do
        BASENAME=${file%%_R1_001.fastq.gz*}

        NUM=`echo $FILE | sed 's/_S.*$//'`

        bowtie2 -x cov_nastya -U $FILE -S ${NUM}exit.sam

done

```

# Answer
It seems you do not have bowtie installed or loaded. Use `module load bowtie2` together when you load `samtools`.
> 3 votes

---
Tags: samtools, bowtie2, slurm
---


title: Obtaining Whole Genetic Sequence

text:
# Question
Title: Obtaining Whole Genetic Sequence
As an end-user of my own data, specifically raw DNA sequence (WGS, whole genome sequencing). How or where do one obtain such DNA data so that a biology hobbyist can perform bioinformatics on this?

What keywords is best describe to the search engines when finding these online WGS websites).

It would be nice to have a USB stick mailed to me after such DNA HomeKit. What file format is most commonly offered? (VCF?)

Am I looking too far ahead into the future?

# Answer
> 2 votes
Are you looking too far ahead into the future? No. This is certainly possible now. The gold standard being Illumina or BGI short read WGS. Long read sequencing can capture some extra data but is very noisy.

How or where do one obtain such DNA data? There are now many companies. You want one the uses NGS with reads at least 150 bp and garantees at least a mean of 30x read depth. In my opinion you also want one that gives you ownership of your data. Ideally they will give you; fastqs (raw reads), a BAM (mapped fastqs) and a VCF (variant file). I good company will ask you if you want to donate your data to science. Be warned these are not defaults. Many (most?) companies won't give you your raw data citing medico legal risk and then they will on sell it to a pharma companies. I think these guys are OK (as they are in Europe with good privacy laws) https://www.dantelabs.com/ but plz do your own research.

Last, these data could include info pertaining to your health and your family's health. Always good to speak to a genetic counsellor to understand these things.

# Answer
> 3 votes
"DNA data" comes in several forms dependent on the technology used to produce them.

Companies like 23andMe are using SNP chips and those are available for an only rather a limited number of species. To be completely honest, I don't even know what would be the raw data they would provide, but I would suspect it could be a vcf file.

If the chosen technology would be DNA sequencing. DNA sequencing is provided by a few dozen of sequencing centres around the world, some of them associated with universities. DNA sequencing generates files that contain individual sequencing reads with assigned quality scores ("how much the sequencer trusts each base"), the typical formats are `fastq` (Illumina) or `.bam` (long read technologies). Usually, you would also need to find a suitable reference in one of the public databases, for instance, the NCBI portal for human reference. Then map your reads and do variant calling. Of course, then it depends on what would you be interested in.

If you don't care for your own DNA and you are more interested in life around you, you can check ncbi for all sorts of genomic data. If you would be interested to understand how the raw sequencing data work, ebi is a good place to look. Also, in every single paper generating/using sequencing data you there should be mention where in public databases you can find them.

Finally, let me warn you. Some of the genomic analyses are quite computationally heavy (not all though). If you don't have a specific question in mind, it would be wise to start with bacterial sequencing.

Hope you will have fun with DNA data!

---
Tags: vcf, file-formats, wgs
---


# Question
Title: How can I use my Myheritage DNA results file for further analysis?
I had my DNA tested by Myheritage and they sent me a csv file with RSID, Chromosome, position and result (which base) with about 700,000 rows. I understand most analyses of DNA use VCF files but is there anything i can do with this csv file i.e. check for genetic health-related genes? I am not a bioinformaticist but I am a scientist and I can use R and python. Ive heard of the Gnomad database but not sure if i can match things to my csv file.

# Answer
A csv file with RSID, chromosome, position and result is enough for what you want to do and these are the core columns of a VCF (which is just a TSV with some headers describing how it was made).

Given you have ~700k rows I suspect that yes there will be genetic health-related SNVs (single nucleotide variant). Disclaimer; these data could include info pertaining to your health and your family's health, I strongly advise that you speak to a genetic counsellor to understand these things.

Gnomad is a population database so generally speaking (plenty of exceptions of course) if there is a SNV in there more than a few times its probably not causing disease. https://www.ncbi.nlm.nih.gov/clinvar/ and https://omim.org/ are example of disease databases.

Here's a pandas example with clinvar (https://ftp.ncbi.nlm.nih.gov/pub/clinvar/vcf\_GRCh38/ you'll need to match your version).

```
import pandas as pd
df = pd.read_csv(
    'clinvar.vcf.gz',
    sep='\t',
    comment = '#',
    header=None,
    names = ['CHROM','POS','ID','REF','ALT','QUAL','FILTER','INFO']
)
df.head()

CHROM   POS ID  REF ALT QUAL    FILTER  INFO
0   1   930188  846933  G   A   .   .   ALLELEID=824438;CLNDISDB=MedGen:CN517202;CLNDN...
1   1   930203  972363  C   T   .   .   ALLELEID=959431;CLNDISDB=MedGen:CN517202;CLNDN...
2   1   930248  789256  G   A   .   .   AF_ESP=0.00347;AF_EXAC=0.00622;AF_TGP=0.00280;...
3   1   930275  969662  T   G   .   .   ALLELEID=959432;CLNDISDB=MedGen:CN517202;CLNDN...
4   1   930336  843786  G   A   .   .   ALLELEID=824439;CLNDISDB=MedGen:CN517202;CLNDN...

```

Then lets say one of your 700k SNVs is at chrom 1, position 930336 and you have an A

```
print(*df.query('(CHROM == 1) & (POS == 930336) & (ALT == "A")')['INFO'].str.split(';'))

```

Gives a list from which you could pull out CLNSIG which here is Uncertain\_significance

```
['ALLELEID=824439', 'CLNDISDB=MedGen:CN517202', 'CLNDN=not_provided', 'CLNHGVS=NC_000001.11:g.930336G>A', 'CLNREVSTAT=criteria_provided,_single_submitter', 'CLNSIG=Uncertain_significance', 'CLNVC=single_nucleotide_variant', 'CLNVCSO=SO:0001483', 'GENEINFO=SAMD11:148398', 'MC=SO:0001583|missense_variant', 'ORIGIN=1']

```

So you could easily parse your 700k SNVs and get their respective clinical info this way. Please think about this carefully before you do it though, there are many implications including life insurance and mental health! I'm of the opinion that this sort of data is useful but MANY people disagree. Further, this advice is general in nature only, I accept no liabilty for what you or anyone else chooses to do with this publicly available data.
> 1 votes

---
Tags: vcf, phylogenetics
---


# Question
Title: How can I run the latest version of samtools on Slurm?
I am trying to run the following script on slurm to extract information from .sam files. While everything is perfect on my local machine with the samtools version 1.11, it does not work on slurm. I have this error:

```
samtools: command not found

```

That means slurm cannot understand the module if I am not mistaken. What would be a solution here?

```
#!/bin/bash

#SBATCH --time=05:00:00
#SBATCH --account=def-myaccount
#SBATCH --mem=200000M
#SBATCH --mail-user=mymail@gmail.com
#SBATCH --mail-type=All

module load nixpkgs/16.09  intel/2018.3
# module load samtools/1.11 # this is what I tried first
module spider samtools/1.11 # this does not help either...

for FILE in *.sam;

do

        NUM=`echo $FILE | sed -r 's/.sam//' | sed -r 's/nastya//'`

#         echo "$FILE" ;

        samtools flagstat -O tsv $FILE > file${NUM}.txt

done

```

# Answer
> 2 votes
I would start with:

```
module avail

```

Then look to see (a) if there is a samtools module avilable and (b) exactly what it's called.

If it's called samtools/1.11 then:

```
module load samtools/1.11

```

Should work, if it doesn't speak to your sys admin.

---
Tags: samtools, slurm
---


title: Multiple box plot

text:
# Question
Title: Multiple box plot
I am using below code to plot multiple box plot but I am still struggling to understand how I can add whisker to box plot and how to change sample order like in this figure SampleType?

```
library(reshape2)
library(ggplot2)
list.files()
morphData <- read.table("Box.plot.chemical.txt", sep="\t", header=TRUE)
morphData
df.m <- melt(morphData, id.var = "SampleType")

pdf("box.plot.pdf", width = 18, height = 12)

p <- ggplot(data = df.m, aes(x="", y=value)) +

theme_bw() +
theme(
      panel.grid.major=element_blank(),
      panel.grid.minor=element_blank(),
      plot.title = element_text(vjust = -8.5,hjust = 0.1),
      line=element_line(size=1),

      axis.text.x = element_text(size = 10, face ="bold", hjust = 0.5, colour = "black"),
      axis.text.y = element_text(size = 12, face ="bold", hjust = 0.5, colour = "black"),
      axis.ticks.x = element_blank()) +

theme(strip.background=element_rect(fill="white")) +
theme(strip.text=element_text(color="black",face="bold", size =14)) +
geom_boxplot(aes(fill=SampleType),outlier.colour = NA)
p + facet_wrap(.~variable, scales="free")
dev.off();

```

update:

I was able to remove the variable name from the bottom by using above updated code but still getting tick mark on x-axis, and also how I can make these label bold?

update:

It will be a great help if you could also guide me how I can add Tukey's posthoc significance letter in above box plot to make this plot more meaningful. I was able to make bold facet wrap text, now please help me with posthoc test letter. I have updated the code again.

update:

Can you please help to change the order of samples; like I want to change the order of SampleType in this box plot and also how we can remove outlier and will it be possible to add whisker so it will look better. I was able to remove outlier by using outlier.colour = NA, Now It just whisker, but I am not able to understand How I can add whisker to box plot? I have updated the code and figure above.

update:

sorry, I wanted to add errorbar and with this command, it is not adding error bar on all box: `stat_boxplot(geom ="errorbar")`

# Answer
> 1 votes
what do you mean by `How I can put the variable name on the top of each box plot like Rice_P at top and at the bottom as well`? Do you want to keep the label at both top and bottom or remove one of them? From my opinion it's not necessary to use facet in this case, since you only have one variable on x-axis in each facet now. You can just put `variable` as x-axis, and color your box as it is now. Boxes with the same `variable` are going to be grouped together.

You can use `axis.ticks.x` in `theme` to remove the ticks on x-axis similar to `axis.text`. To make label bold maybe you can try `fontface="bold"`

Adding testing letter may be harder, you could check out ggpubr or ggprism package.

you can use `as.factor()` to change sample name to factors and plot them in order. Or you can use `scale_x_discrete(limit = <your sample order>)` to force an order. I don't see outlier in plot. Outliers are usually black dots in box plot. And the box plot you show already has whiskers, so I 'm not sure what you want to add.

I'm not sure whether it's a good idea to plot error bars in boxplot because that is not the convention. But probably you can try using `geom_errorbar` in addition to `geom_boxplot` as another layer but it requires more customization.

---
Tags: r, statistics, ggplot2
---


# Question
Title: How to determine a correct promoter sequence? (Bisulfite sequencing)
I'd like to ask you for help with an university assignment. First let me post the entire question:

> Design primers for bisulfite sequencing of the promoter of the ACE2 gene in microbat Myotis lucifugus (based on the annotation of the Broad Institute Myoluc2.0/myoLuc2). The gene is localized on the contig GL429816. As a transcription start, use the start of the ACE2 transcript variant XM\_023753670.1; once you find ACE2 in the UCSC genome browser .You will find it in NCBI RefSeq track

I don't want anyone to solve it for me but I do have a problem with one aspect. I am not clear on how to determine what is a valid promoter for this gene in microbat. I do know it is the upstream sequence however I am not clear on how to tell the number of upstream bases that are the actual promoter. I did find some sources for the actual promoter sequence determined by much smarter people than myself however those are for ACE in humans so I am afraid that it is not applicable in my case.

I'd like someone to point me in the right direction where I could find my answer would be enough.

# Answer
> 1 votes
I think there is no fixed size. Ideally you would have something like histone ChIP-seq or ATAC-seq data to approximate the relevant upstream interval. Since this is not always possible most people go with -500bp. Others might have different suggestions though. Maybe matching known motifs of transcription factors that regulate that promoter to the upstream region can help approximating a better window. But as Sanger sequencing usually gives about 1kb of sequence I'd just take a 500bp to 1kb window and sequence that, fwd and rev primers for sequencing to have double confirmation.

---
Tags: gene, ngs, methylation, homework, epigenetics
---


# Question
Title: Is there a way to do GWAS on phenotype data that is not normally distributed?
Is there a way to do GWAS on phenotype data that is not normally distributed?

For example, if you were measuring a trait as a proportion, and the majority of the data consisted of 0.00 with a long tail of variation. Or, perhaps the data was skewed the other way where the majority of the observations were 100% with a long tail towards zero.

Is there a GWAS method that is appropriate for these types of distributions?

# Answer
Yes, that can be done.

A common approach for these types of analysis is to carry out a transformation on the data first in order for it to have a normal distribution (the most general approach called an inverse normal transformation, among other names), then run the test on the transformed data.

As an example, for data that are positive and zero-skewed, a log transformation may be appropriate.

The transformation approach does require variation in the data values, though. For a trait that is exactly zero for most values, it might be better to treat it as a binary/categorical trait for the purpose of association testing, rather than a quantitative trait. Most association testing packages should allow associations for binary traits to be tested.
> 2 votes

# Answer
as mentioned in previous answers, transformations are frequently used. One commonly used method is quantile normal transformation. Basically you calculate the quantile from the original data and match it to a standard normal distribution. This transforms your data to a standard normal distribution. Even for traits that are naturally normal distributed like height, it is still a good idea to do this transformation as it has been shown to improve power as well.
> 1 votes

# Answer
Sure, you would use the negative binomial distribution when the majority of data was around 0. In other words the data is clustered with the standard deviation far greater than the mean. Most general linear modelling (GLM) will incorporate this approach.
> 0 votes

---
Tags: genome, statistics, gwas
---


title: co-occuring genes in my list

text:
# Question
Title: co-occuring genes in my list
I have a files which has columns representing genes,sample numbers and sample names (55). The genes are displayed by occurrence in the samples (highest to lowest-55 and 2).I wish to calculate the genes which may co-occur in samples. I tried doing this with `awk 'n=x[$1,$2,$3,$4,$5,$6]{print n"\n"$0;} {x[$1,$2,$3,$4,$5,$6]=$0;}'` adding from 1 to 55 columns, but it doesn't work.

An example here:

```
Nfib    55  A1.S1   A2.S2   A3.S3   A4.S4   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.4   IZ_BH_4.5   IZ_BH_4.6   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Lyn 41  A1.S1   A4.S4   A7.S7   A8.S8   A9.S9   A9.S9*  A10.S10 A11.S11 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.9   IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.7   IZ_BH_4.8   IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

Cdk5rap2    52  A1.S1   A2.S2   A3.S3   A6.S6   A7.S7   A8.S8   A9.S9   A9.S9*  A11.S11 A12.S12 IZ_BH_2.1   IZ_BH_2.2   IZ_BH_2.3   IZ_BH_2.4   IZ_BH_2.5   IZ_BH_2.6   IZ_BH_2.7   IZ_BH_2.8   IZ_BH_2.9   IZ_BH_2.10  IZ_BH_2.11  IZ_BH_2.12  IZ_BH_3.1   IZ_BH_3.2   IZ_BH_3.3   IZ_BH_3.4   IZ_BH_3.5   IZ_BH_3.6   IZ_BH_3.7   IZ_BH_3.8   IZ_BH_3.9   IZ_BH_3.10  IZ_BH_3.11  IZ_BH_3.11* IZ_BH_3.12  IZ_BH_4.1   IZ_BH_4.2   IZ_BH_4.3   IZ_BH_4.4   IZ_BH_4.6   IZ_BH_4.8   IZ_BH_4.9   IZ_BH_4.10  IZ_BH_4.11  IZ_BH_5.1   IZ_BH_5.2   IZ_BH_5.3   IZ_BH_5.4   IZ_BH_5.5   IZ_BH_5.6   IZ_BH_5.7   IZ_BH_5.8

```

The first name is of the gene, followed by number of samples and sample names.

for these 3 genes these are the commonly occuring samples:

```
A1.S1, A7.S7, A8.S8, A9.S9*, A9.S9, A11.S11, IZ_BH_2.1, IZ_BH_2.2, IZ_BH_2.3,IZ_BH_2.4, IZ_BH_2.5,IZ_BH_2.6,IZ_BH_2.7, IZ_BH_2.9, IZ_BH_2.11,IZ_BH_2.12, IZ_BH_3.2,IZ_BH_3.3,IZ_BH_3.4, IZ_BH_3.5, IZ_BH_3.6,IZ_BH_3.8,
IZ_BH_3.9,IZ_BH_3.10,IZ_BH_3.11,IZ_BH_4.1,IZ_BH_4.2,IZ_BH_4.4,IZ_BH_4.8,IZ_BH_4.11,IZ_BH_5.1,IZ_BH_5.3,IZ_BH_5.4,IZ_BH_5.5,IZ_BH_5.6,IZ_BH_5.7,IZ_BH_5.8

```

Kindly help.

# Answer
> 1 votes
I think what makes this especially hard is the way your data is formatted. The problem would be way easier if you had a matrix of presence / absence instead. Each row would be a gene, each column a sample, and the values would be 0 (absent) or 1 (present).

Awk becomes a bit inconvenient for this task, so I'll use python instead. To get such a matrix from your file, you could do the following:

```
import numpy as np
import pandas as pd

# Get total number of genes / samples
all_samples = set()
all_genes = set()
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    all_genes.add(line[0])
    all_samples.update(line[2:])
all_samples = list(all_samples)
all_genes = list(all_genes)

# Initialize presence matrix with correct dimensions
pres_mat = np.zeros((len(all_genes), len(all_samples)), dtype=int)

# Map each sample/gene to a row/col
row_map = {gene: i for i, gene in enumerate(all_genes)}
col_map = {sample: i for i, sample in enumerate(all_samples)}

# Fill matrix incrementally
with open('genes.txt', 'r') as genes:
  for line in genes:
    line = line.split('\t')
    curr_gene = line[0]
    row = row_map[curr_gene]
    for sample in line[2:]:
      col = col_map[sample]
      pres_mat[row, col] = 1

# Add names and write to file
df = pd.DataFrame(pres_mat, columns=all_samples)
df.insert(0, 'gene', all_genes)
df.to_csv('gene_presence.tsv, sep='\t')

```

You should get something that looks like:

```
       gene  IZ_BH_5.5  A8.S8  IZ_BH_2.1  A9.S9*  IZ_BH_5.3  A9.S9  ...  IZ_BH_3.11  IZ_BH_4.7  IZ_BH_4.9  IZ_BH_2.8  IZ_BH_3.5  IZ_BH_2.6  IZ_BH_4.3
0       Lyn          1      1          1       1          1      1  ...           1          1          0          0          1          1          1
1      Nfib          1      1          1       1          1      1  ...           1          1          1          1          1          1          0
2  Cdk5rap2          1      1          1       1          1      1  ...           1          0          1          1          1          1          1

```

With this format, it is now easy to compute the most frequent samples. E.g. `df.sum(axis=0)` to get the number of occurence of each sample for these 3 genes. It would also allow you to do different things, such as computing the pairwise similarity between samples.

---
Tags: awk
---


# Question
Title: Differential expression analysis for a subset of transcripts
I want to check the differential expression of a specific class of transcripts (say, long non-coding RNAs) using DESeq2. Now, I know that the normalization step takes into account the total number of reads per sample. If I only input a subset of all the transcripts, the sequencing depth will change.

So my question is: can I simply input the subset of interest? Or should I input all transcripts, and subset the results later?

# Answer
It's safer to import everything. You'll want all the data for normalization and dispersion estimates.
> 3 votes

---
Tags: deseq2, normalization, subset
---


# Question
Title: extract chimeric and multimap reads from bam file
I am trying to extract all chimeric and multi-map reads from either SAM/BAM file. Is there any simple command to do that? Can I use `htslib` for parsing sam/bam files and extracting based on any flags. Any help on how to start will be appreciated.

# Answer
> 1 votes
This is how I do it.

```
samtools view -@5 -f 0x800 -hb /path/sample.bam > sup.bam

```

Then if you want it as a fasta

```
samtools view -F 0x1 -hb sup.bam | samtools fasta -F 0x1 - > sup.fa

```

---
Tags: sam, htslib
---


# Question
Title: how to get access to genotypes and phenotypes used for a GWAS
I'm a master's student working on genomic prediction of complex traits using deep learning. i'm looking for a dataset of human genotypes and phenotypes that has been used for a GWAS. The only thing i can get from the gwas i found is the summary statistics not the raw data.
Does anyone know if it's possible to get such data (for free ) for a research project?
Thanks.

# Answer
> 0 votes
How about:

If you are looking for more comprehensive information, you could try biobanks - e.g. https://www.ukbiobank.ac.uk. These are usually open (or reasonably priced) for research but the approval process can take time.

A side note - some of GWAS catalogues will have a link to the study where the association originated from, e.g. here. You can sometimes find the datasource there.

---
Tags: genomics, gwas, deep-learning
---


# Question
Title: Making PLINK compatible files from VCF file without phenotype information
I have a big VCF file that I need to convert to, preferably, bed/bim/fam files that are readable by plink. Currently using plink 2. I am aware that this version of plink can be used to convert VCFs into binary peds using

`plink --vcf file.vcf.gz --make-bed --out out`

Or something similar. I also know that in this case, plink will automatically fill the phenotype and sex columns in the bed and fam files with zeros. I know there is a `--pheno` flag that allows specifying a txt file for the phenotypes and a `--make-pheno` flag that can be used for the same effect.

How exactly does this work? How can I encode the phenotype and sex information correctly? Can I just use --make-bed and then manually make a phenotype/sex file to use down the line with the --pheno flag? Or would I do something like

`plink --vcf file --pheno phenotype_file --make_bed --out out_file`

directly? In this case, will there not be a problem with the fact that the generated .bed in the first step does not have this information? What about with sex?

Another option I can think of is just edit the generated .fam to include the correct information and use it downstream, but that information will still not be in the .bed file.

I have no family relationships, so the FID fields would just be 0s (I guess I can add the --const-fid flag for that effect when I'm converting, right?).

# Answer
You can use `--pheno` to specify a phenotype file you want to use, the 1st and 2nd columns being FID and IID. These are just used on the fly and cannot be written to the BED files as I understand. `--make-pheno` is used to write binary phenotypes, which is unsuitable if your phenotype is actually continuous. But I still do not quite understand what you want to do here.
> 2 votes

---
Tags: vcf, gwas, plink
---


# Question
Title: multi--platform GEO analysis based in expression of single gene
*This question has also been asked on Biostars*

I have a challenge in analyzing GSE34941( It has 4 platforms) according to the expression of a specific gene in R. First, which is the best way to group my patients according to the expression of a gene to high- and low-expression? Second, which is the best way to remove the batch effect of multi-platform analysis here?

I really appreciate any help you can provide.

# Answer
> -1 votes
Maybe can try the ComBat to remove the batch effect

---
Tags: r, microarray, geo
---


# Question
Title: Nanopore data clarification regarding merging of samples
We had same set of sample ran on nanopore platform twice.

The issue on hand is.

First run was done it was nearly done but the power went off. Second run was done completely with enough sequencing depth.

Now the question I ran the analysis with the second set which was complete that set of result I have. But so far I have't ran the first part.

**Question1** : if I run the first part do i get major differences? in the result

**Question2** : would it be advisable to merge both run since they are now can be called as twi different batches plus would be bring drastic changes in the output ? And in terms of tool for merging I googled which shows poretools.

# Answer
It is entirely unclear which type of analysis you aim to do. You probably should be more specific when asking questions.

That said, given that the run was going perfectly fine until a power issue stopped it, all data produced should be fine and you can indeed combine those runs.

You probably shouldn't bother with poretools. If you are using fastq files to start your analysis you can just cat the basecalled files together. Also after alignment, you can merge the bam files. If your analysis requires the fast5 files then you can probably just put these in the same folder, but you should tell us which analysis you are doing.
> 1 votes

---
Tags: nanopore
---


# Question
Title: Filtering Sequences (entries) by headers ID from a Fasta file database
First of all, sorry if this question has been posted previously (I could not found a solution accessing the previous Q&A).

I have a fasta file as follow:

```
>LNIV02000036.519060.520603 Pseudomonas aeruginosa
GAACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAA...

>LDLH01000045.361.1876 Pseudomonas aeruginosa
AGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATG...

>LNMQ01000087.5825.7394 Enterococcus faecium
TTTTTATGAGAGTTTGATCCTGGCTCAGGACGAACGCTGGCGGCGTGC...

...

```

I would like to filter sequences or extract entries for a specific species, let's say Pseudomonas (I have others in mind...). I tried a few things but I think the main problem is my skills associated with bioinformatics. I would like to create a file content Pseudomonas entries (~9k entries in the current database). Any suggestion or tip?

Thank you for your time.

Best,
Fabricio

# Answer
> 4 votes
You can do this using `seqkit` as follows:

```
seqkit grep -r -n -p '.*Pseudomonas.*' temp.fa

```

To explain a little, `seqkit grep` will allow you to search FASTA/Q files by sequence name or sequence itself. In this instance:

* `-r` tells that the pattern is a regular expression
* `-n` to match by full name instead of just id
* `-p` to specify the regular expression pattern to search

There are a bunch of other tools that come with `seqkit` that can prove indispensable for basic sequence analyses.

# Answer
> 2 votes
Using Biopython:

```
from Bio import SeqIO

for record in SeqIO.parse("input.fa", "fasta"):
    if record.description.split()[-2:] == ["Pseudomonas", "aeruginosa"]:
        print(record.format("fasta"))

```

# Answer
> 0 votes
"describing the entire issue".

I have a database (silva 16S formatted by DADA2 - fasta format) that contains ~300k entries for bacteria (16S rRNA gene). e.g.,

> AB001441.1.1538 Pseudomonas syringae pv. broussonetiae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGG... KR054096.1.1224 Marinococcus halophilus GAAGGTAAAGGCTTCCCCAAGGGCGACGATGGGTAGCCGACTTGA... AB681405.1.1478 Paenibacillus barengoltzii GACGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAGCGGAGT... ...

I would like to retrieve sequences for some genus, e.g., Pseudomonas. After running the seqkit command (mentioned by you above) I ended up with ~9 of Pseudomonas entries in a fasta file (example below).

> AB001441.1.1538 Pseudomonas syringae pv. broussonetiae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA... AB001442.1.1538 Pseudomonas amygdali pv. eriobotryae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA... ...

This pool of 9k sequences of Pseudomonas has several species of Pseudomonas (e.g., P. putida, P. syringe Et Cetera). Since some species of Pseudomonas have multiple sequences in the database (e.g., P. putida has 1.8k entries) I would like to know how many different species of Pseudomonas are there (in the 9k Pseudomonas entries), by names (i.e., in the filtered Psedumonas fasta file). At the end of the day, I would like to say that in this database I have 9k entries of Pseudomonas which represents 'X' number of Pseudomonas species.

What I'm thinking... Retrieve the full ID (AB001442.1.1538 Pseudomonas amygdali pv. eriobotryaeall) or if possible the portion of the ID (headers) that contains the species names (e.g. Psudomonas amygdali pv. eriobotryaeall) and with this, I can list the species names and count them. Does this explanation help to understand what I'm trying to do?

Let me know if you need more info.

Thank you for your time. Fabricio

---
Tags: fasta, database
---


# Question
Title: Are there any databases for gene co-expression or expression pattern clustering?
I am currently working on gene clustering based on co-expression pattern in mouse brain. The problem is I do not have some solid way to test my result. Are there any suggestions for databases indicating groups of genes that co-express (in mouse brain), or even better indicating some similarity metrics between each pair of genes based on expression pattern?

Also, I do want some databases for differentially-expressed genes ;)

PS: By "differentially-expressed" I am referring to a pair of genes whose expression patterns in a single individual differ.

# Answer
> 1 votes
How about:

Other potential databases of interest:

---
Tags: database, public-databases, differential-expression, gene-expression, coexpression
---