text:
How does DeepVariant construct RGB images from DNA sequences?
DeepVariant is a pipeline to call genetic variants from DNA sequencing data.
A major step, before feeding the CNN, is to translate these DNA sequences into images. It's unclear why and how Google constructs the RGB images from the DNA data. Obviously, DNA is a string over an alphabet with the characters: {A, T, C, G}.
It is even hard to understand how the mapping works based on the source code of their unit tests.
In their figure from the paper: A is Red, C is Green, G is blue, and T is Yellow (G+R), but this is still unclear how they construct the 3xNxN image.

EDIT from google's blog:

In this article we will show the six channels in a row, but in
DeepVariant they are encoded as six layers in the third dimension,
giving each tensor a shape of (100, 221, 6) corresponding to (height,
width, channels). The variant in question is always in the center of
each pileup image, here marked with a small line at the top.
Channels are shown in greyscale below in the following order:
Read base: different intensities represent A, C, G, and T.
Base quality: set by the sequencing machine. White is higher quality.
Mapping quality: set by the aligner. White is higher quality.
Strand of alignment: Black is forward; white is reverse.
Read supports variant: White means the read supports the given
alternate allele, grey means it does not.
Base differs from ref: White means the base is different from the
reference, dark grey means the base matches the reference.

See page 3 of https://www.biorxiv.org/content/biorxiv/suppl/2016/12/19/092890.DC2/092890-1.pdf

Actually the paper has made it clear how they did it. You just have to read the supplementary materials closer.

In their figure from the paper: A is Red, C is Green, G is blue, and T
is Yellow (G+R), but this is still unclear how they the 3xNxN image.

In RGB, each dimension is an NxN image. Since you have three dimensions, so it's 3xNxN. The red dimension was used to encode the nucleotide bases. The green dimension was used to encode quality scores. Finally, the blue dimension was used to encode the strand information.

Obviously, DNA is a string over an alphabet with the characters: {A,
T, C, G}.

This is easy. Lot's of ways. You could do a one-hot encoding or what DeepVariant used:
def get_base_color(base):
    base_to_color = {'A': 250, 'G': 180, 'T': 100, 'C': 30}
    return base_to_color.get(base, 0)

Based on the supplementary material as mentioned in a comment by Devon Ryan:

The second phase of DeepVariant encodes the reference and read support
for each candidate variant into an RGB image. The pseudo-code for this
component is shown below; it contains all of the key operations to
build the image, leaving out for clarity error handling, code to deal
with edge cases like when variants occur close to the start or end of
the chromosome, and the implementation of non-essential and/or obvious
functions.

Here is the main function:
WIDTH = 221
HEIGHT = 100

def create_pileup_images(candidate_variants):
    for candidate in candidate_variants:
        for biallelic_variant in split_into_biallelics(candidate):
            start = biallelic_variant.start - (WIDTH - 1) / 2
            end = WIDTH - span_start
            ref_bases = reference.get_bases(start, end)
            image = Image(WIDTH, HEIGHT)
            row_i = fill_reference_pixels(ref_bases, image)
            for read in reads.get_overlapping(start, end):
                if row_i < HEIGHT and is_usable_read(read):
                    add_read(image, read, row_i)
                    row_i += 1
            yield image

def fill_reference_pixels(ref, image):
    for row in range(5):
        for col in range(WIDTH):
            alpha = 0.4
            ref_base = ref[col]
            red = get_base_color(ref_base)
            green = get_quality_color(60)  # The reference is high quality
            blue = get_strand_color(True)  # The reference is on the positive strand
            image[row, col] = make_pixel(red, green, blue, alpha)
    return 5


text:
How to optimize the number of amplicons ordered for a PCR wet experiment with several genomic ROIs?
Imagine we have an experiment for which we would like to minimize the number of purchased primers. Let's look at the layout of a single primer below:

Assume we have a precomputed set of all necessary amplicons and we are trying to reduce the costs and find the minimum set of amplicons that tile all regions of interest (ROIs) with some flanking.
The algorithm seems to work, but I am not sure it is ideal, so I would like to know how to assert it finds the minimal covering set of amplicons and how to improve the runtime complexity.
First, we define the underline representation of a primer/amplicon (I  use the terms interchangeably), and an auxiliary function for maintaining the same colors among different plotting sessions.
from dataclasses import dataclass
import numpy as np

def rand_color():
   return "#" + "%06x" % np.random.randint(0, 0xFFFFFF)

@dataclass(eq=True,frozen=True)
class Primer:
    left_end:  int
    right_end: int
    left_roi: int
    right_roi: int
    color: str
    level: int

We then create a random set of primers:
N = 25; max_coor = 2_000
min_primer_len, max_primer_len = 200, 500
min_required_flank_len, max_required_flank_len = 10, 50
roi_len = 20
primers = []
for i in range(N):
  left_end = np.random.randint(0,max_coor-max_primer_len)
  primer_length = np.random.randint(min_primer_len,max_primer_len)
  required_roi_flank = np.random.randint(min_required_flank_len,max_required_flank_len)
  left_roi = np.random.randint(left_end+required_roi_flank,left_end+primer_length-required_roi_flank-roi_len-1)
  right_roi = left_roi + roi_len
  primers.append(Primer(left_end, left_end+primer_length, left_roi, right_roi, color=rand_color(), level=i))

The algorithm for finding the minimal primer covering set is as follows:
primers.sort(key=lambda p:p.left_end)

def does_a_cover_b(a,b, flank=0):
  return a.left_end< b.left_roi-flank and a.right_end> b.right_roi+flank

p_tmp = primers.copy()

cover = set()
for p in p_tmp:
  add_to_cover_set = True
  for c in cover:
    if does_a_cover_b(c,p):
      add_to_cover_set = False

  if add_to_cover_set:
    cover.add(p)

Now cover is supposed to be a minimal subset of primers that covers all ROIs. To visualize, we use plot_primer_set twice:
import matplotlib.pyplot as plt
from shapely.geometry import Polygon

def plot_primer_set(primer_set):
  plt.figure(figsize=(4,4), dpi=120)
  for c in primer_set:
    p_all = Polygon([(c.left_end,c.level),(c.left_roi,c.level),(c.right_roi,c.level),(c.right_end,c.level)])
    x,y = p_all.exterior.xy
    plt.plot(x,y, markersize=7, marker='o', linewidth=5, color=c.color)
    plt.yticks(np.arange(0, N, 1))
    axes = plt.axes()
    axes.set_xlim([-100, max_coor+100])

plot_primer_set(primers)
plot_primer_set(cover)

Is this the optimal and correct minimal set?
My complexity is $O(N^2)$, but I expected it to be $O(N \log{N})$, what am I missing?

I ended up implementing the following greedy algorithm that should work fairly good for any binary relation even if it is not transitive nor symmetric.
Assume you have the following binary relation (you may plug in any function you want):
def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

Using the auxiliary functions to define the binary relation, draw the graph, and build the graph:
import networkx as nx

def does_a_cover_b(a,b, min_dist=30, max_dist=130):
  return min_dist < b.left_roi - a.left_end < max_dist or min_dist < a.right_end - b.right_roi < max_dist

def draw(G):
  pos = nx.spring_layout(G,scale=2)
  nx.draw(G, pos, with_labels=True, font_weight='bold', node_size=750, width=2)

def build_a_graph(primer_list):
  G = nx.DiGraph()

  for i,p_i in enumerate(primer_list):
    G.add_node(i, data=p_i)

  for i,p_i in enumerate(primer_list):
    for j,p_j in enumerate(primer_list):
      if i == j:
        continue
      if does_a_cover_b(p_i, p_j):
        G.add_edge(i,j)
  return G

we may implement the following algorithm:
p_set = set(primers)
cover_set = set()
covered_primers = set()
while len(covered_primers)<N and len(p_set)>0:
  G = build_a_graph(p_set)
  nodes =sorted(G.out_degree, key=lambda item: item[1], reverse=True)
  node = nodes[0][0]
  absorber = G.nodes[node]['data']
  p_set.remove(absorber)
  cover_set.add(absorber)
  covered_primers.add(absorber)

  for edge in G.out_edges(node):
    primer = G.nodes[edge[1]]['data']
    p_set.remove(primer)
    covered_primers.add(primer)

print(f"Used a cover set of length {len(cover_set)} primers to cover the total primer set with {len(covered_primers)} primers")
plot_primer_set(primers)
plot_primer_set(cover_set)

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to predict with the pre-trained DNABERT model?
I was curious to give DNA BERT a try. This is a BERT (Bidirectional Encoder Representations from Transformers) model that was trained on short (k=3,4,5, or 6) k-mers of DNA. Overall, It is exciting work since it bridges NLP approaches and genetics data.
I span up a pytroch3.6 machine on SageMaker with a GPU and tried to assess the prediction of its pre-trained DNABERT model:
Preparing the data:
!git clone https://github.com/jerryji1993/DNABERT
%cd DNABERT
!python -m pip install --editable .
%cd examples
!python -m pip install -r requirements.txt
!curl -L https://northwestern.app.box.com/shared/static/g8m974tr86h0pvnpymxq84f1yxlhnvbi --output a6-new-12w-0.zip
!unzip a6-new-12w-0.zip
!cp ./sample_data/ft/prom-core/6/dev.tsv /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0

Running prediction based on the pre-trained data:

!python run_finetune.py \
   --model_type=dna \
   --tokenizer_name=dna6 \
   --model_name_or_path=/home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0 \
   --task_name=dnaprom \
   --do_predict \
   --data_dir=/home/ec2-user/SageMaker/DNABERT/examples/sample_data/ft/prom-core/6  \
   --max_seq_length=75 \
   --per_gpu_pred_batch_size=128   \
   --output_dir=/home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/out \
   --predict_dir=/home/ec2-user/SageMaker/DNABERT/examples/result \
   --n_process=4

Output:
11/09/2020 04:35:57 - WARNING - __main__ -   Process rank: -1, device: cuda, n_gpu: 1, distributed training: False, 16-bits training: False
11/09/2020 04:35:57 - INFO - transformers.configuration_utils -   loading configuration file /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/config.json
11/09/2020 04:35:57 - INFO - transformers.configuration_utils -   Model config BertConfig {
  "architectures": [
    "BertForMaskedLM"
  ],
  "attention_probs_dropout_prob": 0.1,
  "bos_token_id": 0,
  "do_sample": false,
  "eos_token_ids": 0,
  "finetuning_task": "dnaprom",
  "hidden_act": "gelu",
  "hidden_dropout_prob": 0.1,
  "hidden_size": 768,
  "id2label": {
    "0": "LABEL_0",
    "1": "LABEL_1"
  },
  "initializer_range": 0.02,
  "intermediate_size": 3072,
  "is_decoder": false,
  "label2id": {
    "LABEL_0": 0,
    "LABEL_1": 1
  },
  "layer_norm_eps": 1e-12,
  "length_penalty": 1.0,
  "max_length": 20,
  "max_position_embeddings": 512,
  "model_type": "bert",
  "num_attention_heads": 12,
  "num_beams": 1,
  "num_hidden_layers": 12,
  "num_labels": 2,
  "num_return_sequences": 1,
  "num_rnn_layer": 1,
  "output_attentions": false,
  "output_hidden_states": false,
  "output_past": true,
  "pad_token_id": 0,
  "pruned_heads": {},
  "repetition_penalty": 1.0,
  "rnn": "lstm",
  "rnn_dropout": 0.0,
  "rnn_hidden": 768,
  "split": 10,
  "temperature": 1.0,
  "top_k": 50,
  "top_p": 1.0,
  "torchscript": false,
  "type_vocab_size": 2,
  "use_bfloat16": false,
  "vocab_size": 4101
}

============================================================
<class 'transformers.tokenization_dna.DNATokenizer'>
11/09/2020 04:35:57 - INFO - filelock -   Lock 140334164036072 acquired on /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/09/2020 04:35:57 - INFO - transformers.file_utils -   https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt not found in cache or force_download set to True, downloading to /home/ec2-user/.cache/torch/transformers/tmpa9hkn2q1
Downloading: 28.7kB [00:00, 12.5MB/s]
11/09/2020 04:35:57 - INFO - transformers.file_utils -   storing https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt in cache at /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - transformers.file_utils -   creating metadata file for /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - filelock -   Lock 140334164036072 released on /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/09/2020 04:35:57 - INFO - transformers.tokenization_utils -   loading file https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt from cache at /home/ec2-user/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/09/2020 04:35:57 - INFO - transformers.modeling_utils -   loading weights file /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/pytorch_model.bin

However, the result directory is empty and I am not sure what did I miss?
!ls -lh /home/ec2-user/SageMaker/DNABERT/examples/result
total 0

!ls -l  /home/ec2-user/SageMaker/DNABERT/examples/6-new-12w-0/out
total 0

EDIT:
I have tried also to run this on another ec2 machine and apparently it gets a segmentation fault, so my conclusion is that it's unstable. I will train a RoBERTAa model or similar by myself:
SageMaker doesn't propagate the segmentation fault to the notebook which is another weird thing about it.
11/12/2020 14:01:18 - INFO - filelock -   Lock 139866003860056 acquired on /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/12/2020 14:01:18 - INFO - transformers.file_utils -   https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt not found in cache or force_download set to True, downloading to /root/.cache/torch/transformers/tmpwpxfkxlm
Downloading: 28.7kB [00:00, 20.2MB/s]
11/12/2020 14:01:18 - INFO - transformers.file_utils -   storing https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt in cache at /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - transformers.file_utils -   creating metadata file for /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - filelock -   Lock 139866003860056 released on /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e.lock
11/12/2020 14:01:18 - INFO - transformers.tokenization_utils -   loading file https://raw.githubusercontent.com/jerryji1993/DNABERT/master/src/transformers/dnabert-config/bert-config-6/vocab.txt from cache at /root/.cache/torch/transformers/ea1474aad40c1c8ed4e1cb7c11345ddda6df27a857fb29e1d4c901d9b900d32d.26f8bd5a32e49c2a8271a46950754a4a767726709b7741c68723bc1db840a87e
11/12/2020 14:01:18 - INFO - transformers.modeling_utils -   loading weights file /home/ec2-user/DNABERT/examples/6-new-12w-0/pytorch_model.bin
Segmentation fault

I wonder if the data science stackexchange would end up being better for this question. I'd want to remove Michaels answer first though, which would be a shame. Could you maybe just repost this over there and provide a link back for context?

I don't why this RNN k-mer Pytorch model is not producing an output. I do know that,
 "rnn_dropout": 0.0,

... is bad news for any ANN (artificial neural network) let alone LSTM RNN, which is a pretty complex model.
Firstly, you have to assign a meaningful drop out rate and 0.5 is okay.
If dropout rate = 0 as an output then the model hasn't worked. A dropout layer is essential to present overfitting. Its a foundational part of deep learning. No-one emperically knows quite why it works, there's a few approximate good guesses, but a dropout is essential.

The other thing to note is what does "n_process" mean? If you are parallelising a GPU across 4-cores, okay its got a lot of layers but even so thats alot of cores. It ain't a CPU and I'd get it working on a single GPU core first.

I see your points. However, this should work presumably out of the box. I have tried to use google colab too, but it was too complicated there compared to SageMaker.

Yep thats Google! Its deep learning is always geeky  :-)  I'm surprised to hear this about colab because it is supposed to be set up for DL, but can vouch 100% thats what  Google developers do. Tensorflow (Google) vs PyTorch (FaceBook) ... is similar, PyTorch is much less geekified.

I was intending to give DNABERT a try as it was published on GitHub. Sadly the developers haven’t replied back to me. Feel free to give it a try and see if it works for you as-is.

I don't know, because I don't know how they are vectorising their k-mers. In my opinion converting DNA to numbers for DL is key to the results. You are using a solid model with a crazy number of layers (12 is alot). Anyway, it was fun to chat.

They use BERT which is a "crazy" model by definition. Here is the paper in case someone wants to see it biorxiv.org/content/10.1101/2020.09.17.301879v1

Thanks @0x90. Sorry I said 'crazy' what I meant was the number of layers was alot. This is a good thing and was (perhaps not now) the basic method for increasing the accuracy of the model in any field of DL. This approach was famously deployed by one big player (was it Google?) to win an image training competition with the other big players. Its just in bioinformatics its the most layers I've seen in any model.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Are there computational tools to extract features of DNA sequences?
I am looking for tools to extract features from short DNA sequences. For example, entropy, complexity, GC-content, etc.
I have found the generateFeatures.py script from the PyFeat repo, but is there love a more widely used source code or a standard way to extract features from sequences biopython or similar?
Also, I think that since there are $2N$ encoded bits in a sequence of $N$ nucleotides, we have at most $2N$ independent features that could be extracted.

In addition, I am curios if there are any transformers models for DNA sequences.

Edit:
In addition, you can use DeepHF's feature utils, which can be accessed and consumed as seen here.

Great question! +1

If you want to do this with biopython, the SeqUtils package could be a solution.


text:
How to compute a p-value for the similarity measurement of two DNA sequences/contigs?
I'm attempting to come up with a statistical model to determine the p-value for the similarity value $ 0 \le s \le 1 $ of two DNA contigs.
The idea is to get a numerical measure to estimate if two contigs are homologous.
I want to account for base composition (and/or other structural features) if possible, but not necessarily because I might want to do the analysis only from the information about the length of the sequences/contigs and largest overlaps. Maybe something inspired from other statistical significance in bioinformatics .
Assume we have two DNA sequences: $a$ and $b$ with the basepair lengths $a_{\text{len}}$, $b_{\text{len}}$. Their overlap $o_{a-b}$ is the number of contiguous bases that appears in both sequences. If we define the similarity $s = \frac{o_{a-b}}{\min{(a_{\text{len}}, b_{\text{len}})}}$ how can we define a p-value for $s$?

Okay, can you please give some more context in your question about a specific use case for this algorithm? It will help to get better answers.

What is basis of your denominator ? How will you calculate ? what does  ,  comma here, stand for? How is it better measure ?

Well, here my thoughts:
The equation $s$ somehow tries to represent the ratio $\frac{equality}{equality+inequality}$. But $s$ returns only a value. And to assing the p-value to a single value... is complicated.
How could you approach this?

Using some kind of Hamming distance for each base pair at the position $n$ ($1 \leqslant n \leqslant len$) between the sequence i ($s_{i}$) and sequence j ($s_{j}$) in order to get a set of binary data:

$$s_{ij} = \begin{Bmatrix} 1\ if\ seq_{i_{n}}=seq_{j_{n}} \\ 0\ if\ seq_{i_{n}}\neq seq_{j_{n}}   \end{Bmatrix}$$
With this, you get a set of binomial data (eg $\{0,1,0,1,0,1,1,1,0,1,1,1,0,...\}$). Now, you should be able to perform some (Binary?) test and get the p-value.
PD: Consider aligning the two sequences if there is interest in discussing base pair homology.

The overlap is what the aligner labeled as the largest similar subsequence

Maybe something along the lines of how blast compute p-value


text:
How to access XL-MS raw data and databases?
I'm looking for databases and/or raw data of protein mass spectrometry experiments.
After reviewing xlinkdb and Proxl, I was curious to learn about other alternatives and options that allow users to run the database locally and browse raw experimental data. NCBI and proteomexchange.org don't seem to have enough data either.


text:
How to use hgvs to project a variant list on a given protein sequence string?
I have the following protein structural variance list:
mutations = ['F594_D600dup', 'D586_E596dup', 'D593_F594insSPEDNEYFYVD', 'E598_F612dup', 'E598_Y599insSGSSDNEYFYVDFREYE', 'E598_Y599insVAYVDFREYE', 'E604_F605insSPRGGNEYFYVDFREYEYDLKWE', 'F594_R595insGTGSSDNEYFYVDF', 'F612_G613insQGFYVDFREYEYDLKWEFPRENLEF', 'L601_K602insNVDFREYEYDL', 'M578_E598dup', 'N587_D600dup', 'R595_E596insDYVDFR', 'R595_L601dup', 'S584_D600dup', 'S584_F605dup', 'S585_F594dup', 'Y597_E598insAGSSDNEYFYVDFREY', 'Y597_E598insDEYFYVDFREY', 'Y597_K602dup', 'Y599_D600insEYEYEYEY', 'Y599_D600insPAPQIMSTSTLISENMNIA', 'Y599_E604dup'])

and I have this protein sequence:
seq="MPALARDGGQLPLLVVFSAMIFGTITNQDLPVIKCVLINHKNNDSSVGKSSSYPMVSESPEDLGCALRPQSSGTVYEAAAVEVDVSASITLQVLVDAPGNISCLWVFKHSSLNCQPHFDLQNRGVVSMVILKMTETQAGEYLLFIQSEATNYTILFTVSIRNTLLYTLRRPYFRKMENQDALVCISESVPEPIVEWVLCDSQGESCKEESPAVVKKEEKVLHELFGTDIRCCARNELGRECTRLFTIDLNQTPQTTLPQLFLKVGEPLWIRCKAVHVNHGFGLTWELENKALEEGNYFEMSTYSTNRTMIRILFAFVSSVARNDTGYYTCSSSKHPSQSALVTIVEKGFINATNSSEDYEIDQYEEFCFSVRFKAYPQIRCTWTFSRKSFPCEQKGLDNGYSISKFCNHKHQPGEYIFHAENDDAQFTKMFTLNIRRKPQVLAEASASQASCFSDGYPLPSWTWKKCSDKSPNCTEEITEGVWNRKANRKVFGQWVSSSTLNMSEAIKGFLVKCCAYNSLGTSCETILLNSPGPFPFIQDNISFYATIGVCLLFIVVLTLLICHKYKKQFRYESQLQMVQVTGSSDNEYFYVDFREYEYDLKWEFPRENLEFGKVLGSGAFGKVMNATAYGISKTGVSIQVAVKMLKEKADSSEREALMSELKMMTQLGSHENIVNLLGACTLSGPIYLIFEYCCYGDLLNYLRSKREKFHRTWTEIFKEHNFSFYPTFQSHPNSSMPGSREVQIHPDSDQISGLHGNSFHSEDEIEYENQKRLEEEEDLNVLTFEDLLCFAYQVAKGMEFLEFKSCVHRDLAARNVLVTHGKVVKICDFGLARDIMSDSNYVVRGNARLPVKWMAPESLFEGIYTIKSDVWSYGILLWEIFSLGVNPYPGIPVDANFYKLIQNGFKMDQPFYATEEIYIIMQSCWAFDSRKRPSFPNLTSFLGCQLADAEEAMYQNVDGRVSECPHTYQNRRPFSREMDLGLLSPQAQVEDS"

I would be happy to work with hgvs in order to project the mutations on the given sequence. However I am not sure how to feed hgvs with an inline sequence (not via a reference sequence)
First, we need to convert all single letters to a three-letter format. So for example, E27del will be converted to Glu27del
str = "Test:p.Glu27del"
hp = hgvs.parser.Parser()
var_str = hp.parse_hgvs_variant(str)
print(var_str)

But I am not sure how to project the mutation if the sequence is simply a string. In the example from their paper they use hg19 as a reference and I want to inject a short sequence instead:
import hgvs.variantmapper
vm = hgvs.assemblymapper.AssemblyMapper(
    hdp, assembly_name='GRCh37', alt_aln_method='splign')
var_g = vm.c_to_g(var_c1)
var_g


text:
How to get the genomic sequences from a blat result?
Assume we have a query.fa file that contains sequences and we run:
blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
the output output.pslx file looks like this:
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block       blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,
20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,

What would be a reasonable way to get the genomic contexts (5bp upsteam and 5bp downstream) for each aligned sequence.
For example, assume that blat found that the seq: AAATTGGGGAAAA aligns to chr2:100-113, so the question is how to get chr2:95-118 easily without reinventing the wheel.

I couldn't make it work with bedtools, because my genome's index file is corrupted, but this should work for others who have successfully used bwa or samtools to index their reference genome:
blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
awk 'NR>5 {print $14 "\t" $16-10"\t" $17+10}' output.pslx > regions.bed
bedtools getfasta -fi /genomes/mm10.fa.gz -bed regions.bed

https://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html

Creat a bed file with chr2 95 118 inside, and then use it with bedtools getfasta to extract your region

@user3479780 I tried that. Add some issue with indexing my reference genome. But this should work too. You are right.

Via BEDOPS convert2bed (psl2bed) and bedops operations:
$ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed

The file answer.bed will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.
This BED file can be run through samtools faidx or similar to get sequence data.
References:

https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range

I couldn't get it to work with varied: pcl2bed, samtools, bedtools commands. I think this biopython should work:
import gzip
from Bio import SeqIO
n=5

exact_matches = {}
blat_qresult = SearchIO.read('blat_output.pslx', 'blat-psl', pslx=True)
for hit in blat_qresult:
    for f in hit.fragments:
        if f.hit_span == f.query_span:
            chrom = hit.id
            if chrom not in exact_matches:
                exact_matches[chrom] = []
            exact_matches[chrom].append(f.hit_range)

with gzip.open('/genomes/mm10.fa.gz', "rt") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        if record.id in exact_matches:
            for subseq in exact_matches[record.id]:
                print(record.seq[subseq[0]-n:subseq[1]+n])


text:
Constructing a phylogenetic tree for classical MHC genes from different organisms
I am interested in constructing a phylogenetic tree containing the classical HLA α1 and α2 domains and potential classical MHC α1 and α2 domains of a different organism. I am using the amino acid sequences for the tree.
However, I am unsure of what bioinformatic tools and software to use and what parameters to fill in.
Lastly any advice on how to read a phylogenetic tree and how to relate the sequence similarity between two sequences by looking at the tree would be very helpful.
My aim is to investigate the classical MHC genes of a non human organism.
Any help and advice would be greatly appreciated!

This is a very unusual situation as a main effort is finding the outgroup where these two domains aren't duplicated.
General situation
In animals a gene tree nearly always follows a species tree, with the exception of gene duplication events, which frequently stem from genome duplication events. Making a tree with lots of genes gets hard to follow and actually becomes pointlessly problematic when you get a branch that does not quite follow the species phylogram —say, you have rabbit, mouse and human genes and the inferred tree, due to poor signal, comes out with mouse+human together as opposed to rabbit+mouse. As a result a nice tree needs to present the species before the genome duplication, often sea lamprey —I wrote a blog post recently were I discuss nice species to pick for a multiple sequence alignment.
In this simple case, one can get the protein sequence of interest and blastp it against human, sea lamprey and fruit fly, see how divergent are the paralogues in human and if the split is more recent that the latter two, if so, the latter will be an outgroup and after getting a few representative species in between, one can make a MSA with Muscle, which contains also a tool to make a neighbour joining tree.
Domain duplication
In your case you have something more complicated in that you have a series of paralogues, each gene encoding a protein (e.g. HLA-A) containing two globular domains, which form the MHC-1 fold, and you want make a tree of this... which means you want find out when the two domains duplicated (i.e. outgroup)? If the domains work together (correct?) and have the same fold (as you alluded to, correct?), this would mean that before the duplication event of the domains, the ancestral protein may have worked as a homodimer.
I am not familiar with the literature, so do not know if this has already been looked at —but I suspect it has.
If not, clicking on the sequence of one of the two domains in the Uniprot entry of one of the paralogues and doing a few iterations of psi-blast against human and a few other distant species (such as zebra fish, sea urchins, which are basal to vertebrates, and fruit fly, which I do not believe have a MHC etc.) most of these will have two hits per protein... but one might not. This is the outgroup!


text:
Download program for ΔG for TM insertion
I am trying to calculate the ΔG of TM insertion for many protein sequences.
Has anyone come across the source code or a programmatic version for the von Heijne ΔG for TM insertion? The web server is at http://dgpred.cbr.su.se/index.php?p=instructions
I can't find a link in the papers or on the site.

Ask the authors otherwise you'll require a "digger", (web spider). I'll share my "digger" code for another delta G website but I have not used this specific site. I think it is common for web servers not to distribute the "backend" for delta G calculations. If you get the standalone please post back. I'm about to start a very similar project (not on TM, but stability of a single antigen) and coding to explore the delta G parameter space is of a lot of interest

@MichaelG. Contacting the authors is the next step. I just want to check that I'm not missing their repo. A digger/scraper isn't ideal, unfortunately. The tool I'm building needs to be a bit of futureproof.

The application is stored in github (ElofssonLab/dgpred) with a link to the website.
To find it I searched for: dgpred download github


text:
Is there any way of using biopython to write Swissprot files?
I have the following mwe for filtering a Swissprot file based on a certain feature, in this case, transmembrane regions.
from Bio import SeqIO
records=[]
for record in SeqIO.parse("Input.txt", "swiss"):
    transmembrane_protein=False
    print record.id
    for i, feature in enumerate(record.features):
        if feature.type == "TRANSMEM":
            transmembrane_protein=True
    if transmembrane_protein==True:
        records.append(record)
SeqIO.write(records, "Output.txt", "swiss")

The script works when SeqIO.write(records, "Output.txt", "swiss") becomes SeqIO.write(records, "Output.txt", "fasta")
However, such a method is not yet supported.
ValueError: Reading format 'swiss' is supported, but not writing

From the docs, I see that writing is not supported for swiss:

Note that while Bio.SeqIO can read all the above file formats, it cannot write to all of them.

Is there any unofficial way of using biopython/python to write swissprot files from parsed swissprot files?

I read that "SeqIO cannot output the file format (e.g. the plain text SwissProt format, "swiss") ". But you could paste the proteins without using Biopython once you know the id of the features you want

It must be in some of the docs somewhere I'm sure, I just haven't seen it explicitly stated, but it is pointed out in the error handling so that's good enough for me to know there is no easy way. I'm dealing with several proteomes. Manual pasting wouldn't be fun ;) I could really do with a programmatic solution.

You could use grep, or awk to select the records and output them in another file using a pipe.

@Llopis That would certainly be one way forward. However I'd like to keep the script in a single python file because it's cross platform and transportable.

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry:
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.
from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())

Perfect solution, worked like a charm.

I don't know how to do this with biopython, but it's simple enough in Perl using Swisskinfe:
#!/usr/bin/env perl
use strict;
use warnings;
use SWISS::Entry;

my $flat=$ARGV[0]||die("Need a flat file\n");

## Change the line termination string so we read an
## entire entry at a time
local $/ = "\n//\n";

my $fh;
## Deal with compressed or uncompressed files
if($flat=~/\.gz$/){
    open($fh,"zcat $flat |")|| die("cannot open flat file $flat : $! \n");
}
else{
    open($fh,"$flat")|| die("cannot open flat file $flat : $!: $@\n");
}
## Parse the file.
entry: while (my $current_entry = <$fh>) {
    ## Read the current entry as an "entry" object
    my $entry = SWISS::Entry->fromText($current_entry);
    ## Iterate over the entry's features
    foreach my $feature ($entry->FTs->elements()) {
        ## If any of the features are 'TRANSMEM', print the entry
        if ($feature->[0] eq 'TRANSMEM') {
            print "$current_entry\n";
            ## Skip to the next entry
            next entry;
        }
    }
}

You need to install the Swissknife module first:

Download Swissknife_1.73.tar.gz from here (SourceForge link).
Install
tar xvzf Swissknife.tar.gz
cd Swissknife_1.73/
perl Makefile.PL
sudo make install

Now, make the script above executable and run with:
foo.pl file.flat

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
How to quickly and robustly convert between mmCIF and PDB?
There is already a question on PDB/CIF to MMTF, however what is a robust way to programmatically go between PDB and CIF files?
For example I can use a python script from this gist that relies on Biopython.
id="4ckh"
wget http://www.ebi.ac.uk/pdbe/static/entry/download/$id-assembly-1.cif.gz
gunzip $id-assembly-1.cif.gz
python cif2pdb.py $id-assembly-1.cif $id.pdb

This process yields a pdb structure as shown below.

The issues are:

Requires python/Biopython to be installed.
The output isn't perfectly clean.
The gist is relatively undeveloped and hasn't got much evidence of testing.

I can't really address your issues certainly so I thought I would comment rather than answer, but I've had good results using Maxit (https://sw-tools.rcsb.org/) which is produced by the RCSB. At least for submitting to the PDB Maxit is endorsed so I would hope it's robust.

Converting a PDB file to some version of an mmCIF file is always possible provided the PDB file doesn't have any "problems". However grouping identical molecules into the same entity is non-trivial so the Biopython and BioJulia mmCIF writers (which I wrote) do not try and do this, simply assigning new entites for each molecule. These files will also be missing other things like the metadata parts of the mmCIF dictionary not involving the atoms.
In Biopython:
from Bio.PDB import PDBParser, MMCIFIO
p = PDBParser()
struc = p.get_structure("", "file.pdb")
io = MMCIFIO()
io.set_structure(struc)
io.save("file.cif")

In BioJulia:
using BioStructures
struc = read("file.pdb", PDB)
writemmcif("file.cif", struc)

Going the other way is, however, not always possible. You can run out of atom serials or chain IDs to use in the output PDB file (this is one reason for the switch to mmCIF files in general). The gist you have linked is a way to get round this.
If this is not a problem for your structure you can invert the above snippets to get a decent PDB file - see the relevant docs. Note that in both cases it will use the auth_seq_id and auth_asym_id fields to get the residue number and chain ID rather than the relabelled label versions. The test suites of Biopython and BioJulia do test these conversions pretty thoroughly.
There is also pdb-tools if you do not want to install Biopython, though both involve a simple pip install so it's not much difference.
For some speed benchmarks of various software see here and here.

hey @jgreener! Thanks for the clarification!
I'm wondering how come MMCIFIO drops metadata and such? Isn't that pretty inconvenient? let's say i need to rename (change the ids of ) certain chains in a structure and nothing else... I'll be losing all the authors' annotation, seq data etc. just by doing io.set_structure(struct), right? Are there any alternatives that you know of?

Yes, you would lose that info. The readers and writers focus on parsing the atomic structural data into a unified representation that can be modified and written back out to any file type. Developer focus has gone there rather than on writing a unified representation for the metadata. My recommendation would either be to concatenate the newly generated file onto the old metadata (beware some may be wrong now) or to use the MMCIF2Dict constructor and modify list items (writing that out preserves metadata, again some metadata may now be wrong).

that's super useful. I have played with MMCIF2Dict before and that'd actually be ideal to preserve and edit the metadata. Any tips for how to write that representation (the dict resultant from MMCIF2Dict) to a .cif file if not through MMCIFIO?

You can write a MMCIF2Dict with MMCIFIO, just use io.set_dict(d) rather than set_structure.


text:
How to calculate Fst from AMOVA
I calculated an AMOVA from a genind object, with one hierarchical factor. In the table I obtain there are SSD values (for my grouping factor,"Error" and total) and sigma2 values (for my grouping factor and "Error")
I have two questions:

What does "error" stand for?
How do I calculate Fst? Which values do I have to use?

A genind object in R is an object which contains allelic information about a set of individuals.
This is the table that I get:
Analysis of Molecular Variance

Call: amova(formula = plantula.dist ~ g, nperm = 100)

           SSD         MSD  df
g     354992.9 354992.9297   1
Error 310814.9    693.7834 448
Total 665807.9   1482.8683 449

Variance components:
       sigma2 P.value
g     1574.69       0
Error  693.78

Variance coefficients:
       a
224.9956

Welcome to Bioinformatics! Could you please clarify what is a geneind object? And could you provide the exact output of the table? (If you provide how you get that table that would probably be also helpful to answer you)

I tried to insert the table, though I can't quite format it right... sorry, I'm still figuring this out!

By Fst are you refering to Fixation index?

Please also indicate the non-standard packages you are using. I’m assuming amova comes from ‹pegas›?

Yes @Llopis I mean fixation index, and yes Konard Rudolph, the amova I am using is from the package pegas.
Sorry for not specifying!

Error is the estimation that can't be fitted by your parameters (in this case g).
The amova function in pegas package does the following:

This function performs a hierarchical analysis of molecular variance as described in Excoffier et al.
  (1992). This implementation accepts any number of hierarchical levels

And as per the definition of Fixation index, you need the variance in the frequency of the allele between different subpopulations, weighted by the sizes of the subpopulations, the variance of the allelic state in the total population. Which is not provided by amova. Thus you cannot compute the fixation index from amova results.

That's what I thought. Thank you for confirming it. And thanks for the explanation for error!

@BeatriceBaldi if this answers your question, you can mark as solved by checking the tick at the left of the answer. This way you let others know that this problem is solved

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
How to troubleshoot access to tabix file via http or https?
I have some BED files that I have compressed with bgzip and indexed with tabix.
I am using version 1.9-392-gcfa32f7-dirty on the server side.
$ /usr/local/bin/bgzip --version
bgzip (htslib) 1.9-392-gcfa32f7-dirty
Copyright (C) 2019 Genome Research Ltd.

The files are served from a subdirectory made available via nginx on port 80 (http) and 443 (https).
The problem is that, for some bgzip-compressed files, I get a correct response to a tabix query if I run it locally, e.g.:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KL.gz" chr1:50669400-50669800
chr1    50669200    50669400    id:253347,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]
chr1    50669400    50669600    id:253348,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]
chr1    50669600    50669800    id:253349,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.1764,15] ]

However, when I run this command from another host (such as my laptop), I get no result:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KL.gz" chr1:50669400-50669800
$

Running echo $? gives a result of 0 (no error). There is nothing in nginx error logs to suggest an issue with http access.
Other bgzip-compressed files work, which further suggests a problem with either the compressed file or its tbi index file, and not with the nginx configuration.
From my laptop, for instance, the following file can be accessed:
$ tabix "http://explore.altius.org/tabix/vB/hg19.15.all.KLs.gz" chr1:50669400-50669800
chr1    50669200    50669400    id:253347,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]
chr1    50669400    50669600    id:253348,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]
chr1    50669600    50669800    id:253349,qcat:[ [0,1], [0,2], [0,3], [0,4], [0,5], [0,6], [0,7], [0,8], [0,9], [0,10], [0,11], [0,12], [0,13], [0,14], [0.23635,15] ]

I've investigated the two test files and they are both sorted lexicographically by chromosome name and numerically by start and stop coordinates. I do not see any odd line endings. These two test files are both virtually identical, except for the data in the fourth column.
Is there anything I can do to further troubleshoot this and figure out why I cannot get a result from other hosts other than the server itself, and only for certain files? (Is there a debug version of htslib applications that would include verbose logging, say?)

Update
In response to Jonathan Moore's comment, I wget the bgzip-ed gz archive and its associated tbi index, and then try to query the file locally (not via http, but directly via the filesystem).
I get the following error message:
$ tabix hg19.15.all.KL.gz chr1:50669400-50669800
[W::hts_idx_load2] The index file is older than the data file: hg19.15.all.KL.gz.tbi

The tbi index file was made immediately after creating the gz file, so I don't understand this error.

Can you wget these particular files from the laptop?  If so, can you then use tabix locally on the laptop to decode the file you just got?

When I do so, I get a strange error. I have updated my question with more details.

That looks to be a fairly common warning message that people get in a few different strange circumstances.  The only thing that comes to mind is whether there is a weird version difference between the tabix format installed on the laptop and that on the server that is coming into play.  Can you get the tabix's to tell you their version numbers?

My other thought is, if you generate a tabix index on the laptop, is it different from the one that you wget from the server?


text:
Tabix queries spanning chromosomes
I can do something like this:
$ tabix ${url} ${range}

where ${range} spans a region on a single chromosome, e.g., chrA:X-Y.
Can a tabix query span multiple chromosomes, e.g., where ${range} looks like chrA:X-chrN:Y (or similar), which returns data between chrA:X and chrN:Y (using either a lexicographical or natural sort order between chrA and chrN)?
Note: I do not want to pass in a BED file of intervals. If it can't be done without passing in a separate file, that's fine. I only want to provide a one-liner to give to someone who can run tabix on their end, to do this type of query directly and remotely.

Looking at the htslib source code, it seems like the answer is no, a single region specifier cannot span multliple chromosomes.
Seems like the closest you can get is to enumerate all the chromosomes as separate regions on the command line. So to get everything between chrA:X and chrN:Y, you would have to write:
tabix inputfile.gz chrA:X chrB chrC ... chrN:1-Y

(If you omit the end coordinate like chrA:X, the region encompasses chrA from the specified coordinate to the end. If only specify the chromosome name without the coordinates, like chrB, the region encompasses the whole named chromosome, and for the last chromosome, you can specify up to the desired point by starting from the beginning: chrN:1-Y.)

Okay, thanks, that worked great. I was able to order chromosomes as desired to write standard output in BEDOPS sort-bed order.


text:
gtf2bed only retains gene information
I want to convert an hg38 GTF file to bed format using gtf2bed:
gtf2bed Homo_sapiens.GRCh38.106.gtf > Homo_sapiens.GRCh38.106.gtf.bed

It outputs a bed file but it only contains the coordinates for the whole genes and has lost the locations of the exons, UTRs, etc.
While I am running gtf2bed, I get thousands of warning like this:
WARNING:root:None   21  33276625    33276753    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33279918 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33279752    33279918    +   three_prime_utr
WARNING:root:WARNING End is after gene end (33283212 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  33283094    33283212    +   three_prime_utr
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   transcript
WARNING:root:WARNING End is after gene end (37101343 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37100814    37101343    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37237744    +   transcript
WARNING:root:WARNING End is after gene end (37221597 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37221419    37221597    +   exon
WARNING:root:WARNING End is after gene end (37237744 > 0)
WARNING:root:WARNING Start is after gene end
WARNING:root:None   21  37237404    37237744    +   exon

How can I create a bed file that also contains the coordinates of the other genomic features?

I cannot reproduce this. Be sure to use the latest version of BEDOPS and to use the correct command:

wget http://ftp.ensembl.org/pub/release-106/gtf/homo_sapiens/Homo_sapiens.GRCh38.106.chr.gtf.gz

gtf2bed < <(zcat Homo_sapiens.GRCh38.106.chr.gtf.gz) > Homo_sapiens.GRCh38.106.chr.bed

$ head Homo_sapiens.GRCh38.106.chr.bed
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  gene    .       gene_id "ENSG00000186827"; gene_version "11"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding";
1       1211339 1214153 ENSG00000186827 .       -       ensembl_havana  transcript      .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214153 ENSG00000186827 .       -       ensembl_havana  exon    .       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; exon_id "ENSE00001832731"; exon_version "2"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1213982 1214127 ENSG00000186827 .       -       ensembl_havana  CDS     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; protein_id "ENSP00000368538"; protein_version "3"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";
1       1214124 1214127 ENSG00000186827 .       -       ensembl_havana  start_codon     0       gene_id "ENSG00000186827"; gene_version "11"; transcript_id "ENST00000379236"; transcript_version "4"; exon_number "1"; gene_name "TNFRSF4"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "TNFRSF4-201"; transcript_source "ensembl_havana"; transcript_biotype "protein_coding"; tag "CCDS"; ccds_id "CCDS11"; tag "basic"; transcript_support_level "1 (assigned to previous version 3)";

Finishes without any of the mentioned warnings.

Thanks! I did not know I should have run it as follows: gtf2bed < Homo_sapiens.GRCh38.106.gtf > Homo_sapiens.GRCh38.106.gtf.bed


text:
Run kallisto iteratively across many samples
I am on a Mac using UNIX. I am trying to use the kallisto quant command on all files in a directory (instead of manually entering them). Because I'm running the analysis against the same index file, I first enter the following:
./kallisto index -i --index --make-unique index.fa

This successfully creates an index file. Then, I tried this:
for file in *.fasta; do kallisto quant --single -l 200 -s 0.1 -o $file-aligned
"$file"; done

While the $file-aligned folders are created, they're empty. I get an error that states that the index file is missing. So I assumed I would need to specify the index file in the command line. I then tried this:
for file in *.fasta; do kallisto quant -i index.fa --single -l 200 -s 0.1 -o
$file-aligned  "$file"; done

I get an error that says "Error: incompatible indices. Found version 0, expected version 10.
I have now also tried the following:
for file in *.fasta; do ./kallisto quant -i index --single - 200 -s 0.01 -o
"${file}-aligned" "${file}";

but I get an error saying the index file is not found.

I am now running
or file in *.fasta; do ./kallisto quant -i index --single -l 200 -s 0.1 -o $file-aligned

"$file"; done
This sort of works. The index file error messages are gone, but it runs quant on the same file 10 times instead of each file once. If I run echo $file it spits back the file that my command is running repeatedly.
If I close terminal, start a new session and run the same command, I get
Error: file not found ""

Why are you using fasta files for reads? Kallisto expects FASTQ files as input.

Does kallisto use quality scores of reads? I had an impression that it does not, but I don't have a reference to backup it....

I down-voted because this post does not pose a question.

Running the jobs in parallel
Iteratively is one solution, but since you're not using bootstrapping, which can be run multithread (and is not needed if you're taking the popular tximport route) you can do better using GNU parallel, this will enable you to run as many jobs as you have execution threads on your system simultaneously.  Note that Kalliso only runs multithreaded with -t if you're using bootstrapping.

You'll need to make the index first:
kallisto index --make-unique -I index.fa.idx index.fa

Next lets define the input FASTQ files and the sample IDs you want to use using a tab-delimited input file like so:
sample_1    sample_1_R1.fastq.gz
sample_2    sample_2_R1.fastq.gz
sample_3    sample_3_R1.fastq.gz

Here the first column is the sample ID which will be used for output and the second separated by a tab character is the input file.
You can produce this using a text editor but as a short-cut you can feed the list of fastq.gz you have into a list via:
ls -1 *.fastq.gz > editme

Next we create a list of jobs to run on our system using a utility script Make_job_list_kallisto.sh which consumes our tab delimited input file list above.  I've used your -l and -s parameters below, note that with paired end data normally these are estimated for you.
#!/bin/bash -eu

[ $# -ne 2 ] && { echo -en "\n*** This script generates jobs for GNU parallel. *** \n\n Error Nothing to do, usage: < input tab delimited list > < output run list file >\n\n" ; exit 1; }
set -o pipefail

# Get command-line args
INPUT_LIST=$1
OUTPUT=$2

# Set counter
COUNT=1
END=$(wc -l $INPUT_LIST | awk '{print $1}')

echo " "
echo " * Input file is: $INPUT_LIST"
echo " * Number of runs: $END"
echo " * Output job list for GNU parallel saved to: $OUTPUT"
echo " "

# Main bit of command-line for job
CMD="kallisto quant -i index.fa.idx --single -l 200 -s 0.1"

# Main Loop
[ -e $OUTPUT ] && rm $OUTPUT
while [ $COUNT -le $END ];
do
    LINE=( $(awk "NR==$COUNT" $INPUT_LIST) )
    # Make file list
    echo "Working on $COUNT of $END Sample ID: ${LINE[0]}, Files ${LINE[@]:1}"
    echo "$CMD -o ${LINE[0]} ${LINE[@]:1}" >> $OUTPUT
    ((COUNT++))
done

We use the script like so:
 Make_job_list_kallisto.sh input_file_list.txt job_list.txt

Now we have a job list for GNU parallel we can run our jobs in parallel like so:
parallel --progress --jobs 4 --joblog kallisto_joblog.txt < job_list.txt

This will run 4 of the kallisto jobs simultaneously (Increase this if you have more threads - don't mind consuming all the IO/CPU threads on your system) producing a nice progress report of those completed, additionally a job log will be written to the file kallisto_joblog.txt.

For paired end data I already have these scripts + similar instructions on GitHub which people might find useful.

@ZincFingers, so -I upper case for specifying the index in my run generation scrip should in fact be lower case -i.  I've tested this and it works for me.  This might be your issue or you might have some other error.  If it's still not working try pasting in a few lines of your jobs list here.

The error I was getting was my fault, I fixed it. Now, when I run the parallel command, I get an error from kallisto saying
Error: missing read files
Error: file srr3884878.gz exists and is not a directory
It says this for all 17 of my files

I suspect you still have an issue with you input file to Make_job_list_kallisto.sh or the jobs it's generating, without seeing the input sample ID file list or the job list it's very difficult to debug.

Check the diagnostic output of Make_job_list_kallisto.sh is sane, I get:
* Input file is: test_for_BS.txt
* Number of runs: 4
* Output for GNU parallel saved to: test_job_list.txt
Working on 1 of 4 Sample ID: test1, Files ATCACG_L004_R1_001.fastq.gz
Working on 2 of 4 Sample ID: test2, Files GTCCGC_L004_R1_001.trimmed.fastq.gz
etc.

I think you might have forgotten to give an output dir or sample ID in the first column of your tab delimited input list you're giving to Make_job_list_kallisto.sh because it looks like kallisto is missing the variable which would be supplied to -o and treating the input fastq.gz file as the output dir, it's then not finding the gz.fastq file as it's trying to use that as the output dir.

You were right! I fixed the input file and it seems to have resolved the issue. Many thanks for your fantastic solution and bearing with me while I trudged through this.

Glad you got it working!

For generating the index, you should be using one of -i or --index, but not both:
./kallisto index --make-unique -i index.fa.idx index.fa

You need to point kallisto to the generated index file (usually has a .idx extension), rather than the fasta file:
for file in *.fasta;
  do kallisto quant -i index.fa.idx --single -l 200 -s 0.1 -o "${file}-aligned" "${file}";
done

Why run it sequentially? Shouldn't this be parallelized, i.e. multiprocessed

You're asking the wrong person; try adding that as a comment to the question. Kallisto quant allows multi-threaded operation with -t <threads>.

I get the same error when I point to the index file
Error: incompatible indices. Found version 0, expected version 10

maybe it's the index command that needs a bit of attention... I'll update the commands to include index generation

I managed to get this to work after fiddling with the index generation:
' for f in *.fasta' do kallisto quant -i index --single -l 200 -s 0.1 -o "${file}-aligned" "${file)"; '
But it appears that it's repeatedly analyzing one of the 10 .fasta files I placed in the directory (also I'm not sure how to get that line to appear as a code)

@gringer - just an additional comment on multithreading: as of version 0.43.1, the -t <threads> parameter doesn't work should you wish to look at the alignments (which kallisto quant can output upon supplying --pseudobam. I usually get around it by running several individual kallisto quant commands concurrently.

@gringer Using -t is fine but it’s actually almost certainly much more efficient to process the files in parallel rather than in a sequential loop.

-t is not going to help, it's only useful if you've enabled bootstrapping, otherwise you can set it as high as you want but kallisto will still run single core.  Multi-threaded execution is used only for concurrently running the bootstrapping.

If you're using tximport for downstream analysis in DESeq2 rather than Sleuth the bootstrapping won't be used as it's not taking into account quantification uncertainty via the bootstrapping procedure, so you need to consider if this is needed.


text:
Spliced vs. unspliced ratios for transcripts in RNA-seq data
Is there a computational tool for measuring what percentage of RNA is spliced in an RNAseq experiment?
I'm not particularly interested in complicated analyses that give ratios for all possible alternative splicing variations. I'd rather have a binary classification. For example: (RNA1 = 35% spliced, 65% unspliced)
Edit: (Additional information)
We're very early in the exploratory stage, so a specific hypothesis is still forming. But we're interested in the sequence content of lncRNAs, many of which are poorly spliced. So we want to get a sense for how prevalent unspliced version of these transcript are.
Because this is currently exploratory, estimations and noise are fairly acceptable. A very basic approach might be to calculate the ratio of reads that map to exon-to-exon junctions, verses those that map from exon-to-intro junctions.
A more advanced approach I'm also considering is Salmon.

Did you ever settle on something? I have the same question as you.

We ended up using aligning with STAR then quantifying with salmon. We gave salmon all known splice variants, including the unspliced annotation. Then we summed all spliced isoforms to get "spliced vs. unspliced".
This pipeline was easy, accurate, and we were really satisfied with the result.

I don't know about tools, but i've used the following python code to calculate the ratio of reads that overlap the 5' or 3' ends of introns or that are spliced. We sum these across all introns in a gene set (we actaully use this for iCLIP analysis to see if RNA binding proteins bind pre-mRNA or spliced RNA).
import pysam
from collections import Counter
from CGAT import GTF, IOTools

def calculateSplicingIndex(bamfile, gtffile, outfile):

    bamfile = pysam.AlignmentFile(bamfile)

    counts = Counter()

    for transcript in GTF.transcript_iterator(
            GTF.iterator(IOTools.openFile(gtffile))):

        introns = GTF.toIntronIntervals(transcript)

        for intron in introns:
            reads = bamfile.fetch(
                reference=transcript[0].contig,
                start=intron[0], end=intron[1])

            for read in reads:
                if 'N' in read.cigarstring:
                    blocks = read.get_blocks()
                    starts, ends = zip(*blocks)
                    if intron[0] in ends and intron[1] in starts:
                        counts["Exon_Exon"] += 1
                    else:
                        counts["spliced_uncounted"] += 1
                elif (read.reference_start <= intron[0] - 3
                      and read.reference_end >= intron[0] + 3):
                    if transcript[0].strand == "+":
                        counts["Exon_Intron"] += 1
                    else:
                        counts["Intron_Exon"] += 1
                elif (read.reference_start <= intron[1] - 3
                      and read.reference_end >= intron[1] + 3):
                    if transcript[0].strand == "+":
                        counts["Intron_Exon"] += 1
                    else:
                        counts["Exon_Intron"] += 1
                else:
                    counts["unspliced_uncounted"] += 1

    header = ["Exon_Exon",
              "Exon_Intron",
              "Intron_Exon",
              "spliced_uncounted",
              "unspliced_uncounted"]

    with IOTools.openFile(outfile, "w") as outf:

        outf.write("\t".join(header)+"\n")
        outf.write("\t".join(map(str, [counts[col] for col in header]))
                   + "\n")

Unfortunately this uses a bunch of libraries you may or may not have, including CGAT (for the GTF parser and IOTools package) and pysam.
Once you've got these statistics you can calculate the "splicing index" as
the log2 ratio of 2 times the number of spliced reads divided by the number of reads overlapping the 3' and 5' ends of introns.

Take this first comment with a grain of salt, since this isn't an area I've worked in much, but: is binary classification possible? If a gene has 3 introns, and 2 are spliced out but 1 is retained, is this "spliced" or "unspliced". My first impression is that an analysis would be a bit more nuanced than a binary classification.
That said, I'm aware of a tool called Keep Me Around for intron retention analysis (preprint, code). I've never used this software, but it looks like it's actively maintained and was created by a research group with a lot of experience and influence in this area, for whatever that's worth. :-)
UPDATE: It looks like the original post was modified since I submitted my answer, but I'll leave the text about binary classification around for posterity. :-)


text:
Production database system for storing vcf records
I’m standing up a production bioinformatics pipeline that’ll be producing standard VCF files. I’d like to store this data in an SQL database that ingest the VCFs and allow me and my team to easily query and visualize many VCFs. Because multiple people will have access to the data, roles/permissions need to be built in. I don’t want to create this system from scratch, but it needs to be both robust and intuitive. Are there any products like this available today?
[Edit]: The answers provided to this question are on the right track, but only provide a small slice of the functionality I'm looking for. For starters, I'm looking for a hosted/managed solution; something potentially like Basespace, but more user-friendly and not quite so bloated.
Additionally, the linked solutions are mostly about getting vcfs into a db, but don't make much headway into providing exploration/visualization tooling.

May be something useful here https://bioinformatics.stackexchange.com/questions/15248/how-to-import-large-bed-gff-vcf-paf-sam-files-into-an-sql-database

Hi @JessimeKirk IMO it is better to ask a question focussed on a single specific issue in most cases. Maybe summarise the answer kindly provided above and focus on a Basespace type solution?

If you are open to non-free solutions, check out VarSome Clinical. (disclaimer: I work for the company behind this, so I am obviously not impartial, but it really is good, honest!).

Like terdon says, you may have to go for commercial solutions if you want layman storage/visualization for VCF files. You're looking to offload bioinformatics work and no free tool can effectively replace something a human being needs to do


text:
UCSC packages/tools ucsc-netchainsubset and ucsc-fatotwobit not working on GNU/Linux
I have installed some ucsc tools on my Mac computer and it worked fine:
Installation
conda create --name ucsc-netchainsubset377 --channel bioconda --yes ucsc-netchainsubset=377

Conda environment activation
source activate ucsc-netchainsubset377

Checking the tool's help
    netchainsubset
netChainSubset - Create chain file with subset of chains that appear in the net
usage:
   netChainSubset in.net in.chain out.chain
options:
   -gapOut=gap.tab - Output gap sizes to file
   -type=XXX - Restrict output to particular type in net file
   -splitOnInsert - Split chain when get an insertion of another chain
   -wholeChains - Write entire chain references by net, don't split
    when a high-level net is encoundered.  This is useful when nets
    have been filtered.
   -skipMissing - skip chains that are not found instead of generating
    an error.  Useful if chains have been filtered.

The same commands on GNU/Linux produced the following error:
netchainsubset
netchainsubset: command not found

even though the tool seems to be installed successfully:
Downloading and Extracting Packages
ucsc-netchainsubset- | 310 KB    | ############################################################################ | 100%
mysql-connector-c-6. | 4.4 MB    | ############################################################################ | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

The same happened with ucsc-fatotwobit and ucsc-liftup and I would not be surprised if this behaviour affected the whole suite of ucsc programmes installed through conda.
Anyone can suggest a solution?
Update 5 minutes after posting
Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.
ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt

Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.
Cheers terdon.
Cheers.

Please [edit] your question and tell us how you installed the tools. It sounds like they're simply not in your $PATH, and you need to either add the directory where you installed them to your $PATH or just call the program using its full path (/path/to/ucsc/netchainsubset).

Short answer: respect the case of the tool name
Longer story below
Thanks to terdon, I printed the PATH when the conda environment was active and listed what was in the folder where conda installed the UCSC tool.
ls /mnt/home2/miska/cr517/anaconda3/envs/ucsc-netchainsubset377/bin
c_rehash         libpng-config      mysql_config    openssl  pngfix
libpng16-config  my_print_defaults  netChainSubset  perror   png-fix-itxt
Then I typed netChainSubset and got the help page, suggesting the tool is installed and working. Interestingly, on Mac, the tool is also called netChainSubset but it seems it can be called with netchainsubsetor NETCHAINSUBset for that matter. The case doesn't seem to matter on mac. To make my code more portable, I am going to follow the case in the bin of where the tool is installed.

Yes, the default file system macOS uses, HFS+, is case insensitive. One of the very few file systems that are. That's why the case didn't matter on your mac. However, most file systems, including those used by Linux, are case sensitive and A is not the same as a.


text:
tRNAscan-SE error: FATAL: Unable to find /usr/local/bin/cmsearch executable
I have downloaded tRNAscan-SE from here. After decompressing and untaring the file, I installed it using:
./configure
make
make install

When I type tRNAscan-SE --help, I get the help page:
tRNAscan-SE 2.0 (December 2017)

FATAL: No sequence file(s) specified.

Usage: tRNAscan-SE [-options] <FASTA file(s)>

  Scan a sequence file for tRNAs
          -- default: use Infernal & tRNA covariance models
             with eukaryotic sequences
             (use -B, -A, -M, -O or -G to scan other types of sequences)

Basic Options
  -E         : search for eukaryotic tRNAs (default)
  -B         : search for bacterial tRNAs
  -A         : search for archaeal tRNAs
  -M <model> : search for mitochondrial tRNAs
                 options: mammal, vert
  -O         : search for other organellar tRNAs
  -G         : use general tRNA model (cytoslic tRNAs from all 3 domains included)
  -L         : search using the legacy method (tRNAscan, EufindtRNA, and COVE)
                 use with -E, -B, -A, -O, or -G
  -I         : search using Infernal (default)
                 use with -E, -B, -A, -O, or -G
  -o <file>  : save final results in <file>
  -f <file>  : save tRNA secondary structures to <file>
  -m <file>  : save statistics summary for run in <file>
               (speed, # tRNAs found in each part of search, etc)
  -H         : show both primary and secondary structure components to
               covariance model bit scores
  -q         : quiet mode (credits & run option selections suppressed)

  -h         : print full list (long) of available options

I then tried to use it with tRNAscan-SE c_elegans.PRJNA13758.WS263.genomic.fa but it gave the following error:
FATAL: Unable to find /usr/local/bin/cmsearch executable
After googling cmsearch, it seems that is is part of the Infernal package. So, I installed infernal with conda using conda install infernal.
I pulled the help from cmsearch successfully:
cmsearch -h
# cmsearch :: search CM(s) against a sequence database
# INFERNAL 1.1.2 (July 2016)
# Copyright (C) 2016 Howard Hughes Medical Institute.
# Freely distributed under a BSD open source license.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Usage: cmsearch [options] <cmfile> <seqdb>

Basic options:
  -h        : show brief help on version and usage
  -g        : configure CM for glocal alignment [default: local]
  -Z <x>    : set search space size in *Mb* to <x> for E-value calculations  (x>0)
  --devhelp : show list of otherwise hidden developer/expert options

Options directing output:
  -o <f>       : direct output to file <f>, not stdout
  -A <f>       : save multiple alignment of all significant hits to file <s>
  --tblout <f> : save parseable table of hits to file <s>
  --acc        : prefer accessions over names in output
  --noali      : don't output alignments, so output is smaller
  --notextw    : unlimit ASCII text output line width
  --textw <n>  : set max width of ASCII text output lines  [120]  (n>=120)
  --verbose    : report extra information; mainly useful for debugging

Options controlling reporting thresholds:
  -E <x> : report sequences <= this E-value threshold in output  [10.0]  (x>0)
  -T <x> : report sequences >= this score threshold in output

Options controlling inclusion (significance) thresholds:
  --incE <x> : consider sequences <= this E-value threshold as significant  [0.01]
  --incT <x> : consider sequences >= this score threshold as significant

Options controlling model-specific reporting thresholds:
  --cut_ga : use CM's GA gathering cutoffs as reporting thresholds
  --cut_nc : use CM's NC noise cutoffs as reporting thresholds
  --cut_tc : use CM's TC trusted cutoffs as reporting thresholds

Options controlling acceleration heuristics*:
  --max      : turn all heuristic filters off (slow)
  --nohmm    : skip all HMM filter stages, use only CM (slow)
  --mid      : skip first two HMM filter stages (SSV & Vit)
  --default  : default: run search space size-dependent pipeline  [default]
  --rfam     : set heuristic filters at Rfam-level (fast)
  --hmmonly  : use HMM only, don't use a CM at all
  --FZ <x>   : set filters to defaults used for a search space of size <x> Mb
  --Fmid <x> : with --mid, set P-value threshold for HMM stages to <x>  [0.02]

Other options*:
  --notrunc     : do not allow truncated hits at sequence termini
  --anytrunc    : allow full and truncated hits anywhere within sequences
  --nonull3     : turn off the NULL3 post hoc additional null model
  --mxsize <x>  : set max allowed alnment mx size to <x> Mb [df: autodetermined]
  --smxsize <x> : set max allowed size of search DP matrices to <x> Mb  [128.]
  --cyk         : use scanning CM CYK algorithm, not Inside in final stage
  --acyk        : align hits with CYK, not optimal accuracy
  --wcx <x>     : set W (expected max hit len) as <x> * cm->clen (model len)
  --toponly     : only search the top strand
  --bottomonly  : only search the bottom strand
  --tformat <s> : assert target <seqdb> is in format <s>: no autodetection
  --cpu <n>     : number of parallel CPU workers to use for multithreads

*Use --devhelp to show additional expert options.

I tried the tRNAscan-SE command again but it failed with the same error message.
So, tRNAscan-SE is looking for cmsearch in the wrong place. How can I make tRNAscan-SE use cmsearch successfully?
Cheers.
Update
Forgot to tell you I am working on Mac OS X and that cmsearch is installed here:
/Users/user/anaconda/envs/python3env/bin/cmsearch
and my python3env anaconda environment is active.

Please [edit] your question and tell us i) what operating system you are using and ii) assuming it is some sort of *nix, the output of which cmsearch

An easy way to build tRNAscan and their dependencies "all Infernal packages" without confusion is to run
conda install -c bioconda trnascan-se

You may try to set the correct path through the infernal_dir variable in tRNAscan-SE.conf.src. You might need to do this before the installation process.
Original file starts with:
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: {bin_dir}

Try the following:
# tRNAscan-SE 2.0
# Configuration File

# default paths
bin_dir: @bindir@
lib_dir: @libdir@/tRNAscan-SE
infernal_dir: "/Users/user/anaconda/envs/python3env/bin"

I've done what you said and it still works but how do I know it's not because my symbolic links are still active? Do you know how I can erase them?

@charlesdarwin You can remove them by  (carefully) performingsudo rm /usr/local/bin/cmsearch and sudo rm /usr/local/bin/cmscan.

tRNAscan-SE 2.0 requires separate installation of Infernal 1.1.2.
The source code and prebuilt binaries can be obtained at
http://eddylab.org/infernal/.
The installation directory of Infernal should be the same as the one
where tRNAscan-SE 2.0 is installed.
When install infernal, run the configure command with the option --prefix <your tRNAscan-SE-2.0 install path>.
If you didn't add option --prefix when installing tRNAscan-SE-2.0, then you should install infernal without option --prefix.

A much cleaner solution is given in @bli's answer.

The error message is telling you that tRNAScan-SE is looking for the cmsearch executable at /usr/local/bin/cmsearch. However, you installed it though conda which apparently placed it at /Users/user/anaconda/envs/python3env/bin/cmsearch.
A simple, if not very elegant solution, would be to make a link to it:
sudo ln -s /Users/user/anaconda/envs/python3env/bin/cmsearch /usr/local/bin/

had to do the same with cmscan and now the code is running.

what does this soft link do?

@charlesdarwin This enables the program to access the program from the place it looks at by default. There might be a cleaner way to do: Changing the default (see my untested answer).

@charlesdarwin as bli said, that just creates a "fake file" at /usr/local/bin/ which points to /Users/user/anaconda/envs/python3env/bin/cmsearch, so that when something tries to open /usr/local/bin/cmsearch, it will actually run /Users/user/anaconda/envs/python3env/bin/cmsearch. But bli's approach of changing the configuration file is a much cleaner solution. Please use that instead.


text:
Find a map/correspondence between two versions of a genome
I am working with two versions of the C. elegans genome. I am finding interesting regions (specifically, tRNA genes) in version 1 and then I would like to know if version 2 also has a tRNA gene in that region. For this, I need to know the correspondence between coordinates in genome 1 and coordinates in genome 2. I think a whole-genome alignment programme might help me, e. g. nucmer or mauve. I have run nucmer on the two genomes and got a delta file but I don't know how to interpret it. I have run show-coords on the delta file and obtained a .coords file. It seems more human-friendly. The first few lines look like this:
/Users/user/Documents/project/output/genome/celegans/ref/wormbase/seq/c_elegans.PRJNA13758.WS245.genomic.fa /Users/user/Documents/project/output/genome/celegans/strain/seq/bristol/scaffolds.fa
NUCMER

    [S1]     [E1]  |     [S2]     [E2]  |  [LEN 1]  [LEN 2]  |  [% IDY]  |  [LEN R]  [LEN Q]  |  [COV R]  [COV Q]  | [TAGS]
===============================================================================================================================
       1      432  |      498      929  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I
       1      432  |     1002     1433  |      432      432  |   100.00  | 15072434 15371751  |     0.00     0.00  | I  I

Does this mean that chromosome I starting at 1 and ending at 432 of version 1 corresponds to chromosome I starting at 498 and ending at 929 and that there is 100% identity between these two 432-long sequences? In this case, it's easy to find a correspondence between bases in version 1 and bases in version 2. However, some lines do not align sequences of the same length:
   1   237139  |     4022   241181  |   237139   237160  |    99.99  | 15072434 15371751  |     1.57     1.54  | I  I

Also, in this sequence position 1 of chromosome I is aligned again? So is it aligning to more than one position in version 2 of the genome? Also one sequence is 237139 bases long and the other 237160 bases long. How do I know the internal structure?
Can anybody guide me, please?
I've also aligned the two genomes using Mauve:
progressiveMauve c_elegans.PRJNA13758.WS245.genomic.fa scaffolds.fa --output scaffolds.mauve > scaffolds.txt

I've tried to import the .mauve file into Python using Biopython but I got an error:
alignment = SeqIO.read('scaffolds.mauve', "mauve")

Traceback (most recent call last):
  File "code/mauve.py", line 4, in <module>
    alignment = SeqIO.read('scaffolds.mauve', "mauve")
  File "/Users/user/anaconda/lib/python3.6/site-packages/Bio/SeqIO/__init__.py", line 712, in read
    raise ValueError("More than one record found in handle")
ValueError: More than one record found in handle

Note that the alignment opened in the Mauve GUI looks fine. Unfortunately, the picture won't load here.
I am working on Mac OSX.
I would appreciate your help in understanding the delta file, resolving the Python error or any other support in finding a function between coordinates in version 1 of the genome and coordinates in version 2 of the genome.

SeqIO.read is for a single record, try SeqIO.parse

NCBI has  REMAP page/tool where you can feed it a file with coordinates of sequence features from one build version, and ask for those features mapped to a different build.  I have only used it to go forward, not backword.  Also, each version of the C. elegans genome used to have a liftover directory at the Sanger ftp site that showed the version changes, as I recall. Check out the WormBase site.

I don't see where I can feed my own file to the Remap service of NCBI: https://www.ncbi.nlm.nih.gov/genome/tools/remap. Reply from NCBI to my request: "The remap service relies on assembly to assembly alignment data to do this remaping/conversion. These datasets came from NCBI's eukaryotic genome assembly workflow. For outside genomic assembly that are not deposited to NCBI, we have no access, so we will not be able to provide mapping between that and those available from the NCBI collection."

The file that you need is a bed file, or a GFF3 file, or similar, that contains the start and end coordinates for all of your tRNA gene regions from genome Version 1.  You select and upload that file and tell REMAP the version/build using one drop down menu. Then in another drop down menu you select the genome version/build for Version 2. Click on the button and in a little while you will have some files to download showing the coordinates for those gene regions in Version 1, and how they map to Version 2.

@mdperry the versions I am using are personal. I have reassembled the genomes using my own data so I cannot choose a version available on NCBI because I have not deposited my versions yet.

Hmmmm, okay; that was not clear to me from your original post.  I see that one of the files you mention has ‘WS245’ in the name and I assumed that was a specific official build of the C. elegans genome, and I erroneously concluded that the other ‘version’ (you mention version 1 and version 2) was a different specific official genome build, and my comments were for that particular situation.  If you are using your own, homemade versions wont you eventually have to map the genes of interest back to one of the official genome builds (i.e., prior to publication)?

That's a good question. For the time being, I use my own build because I think it's better than WormBase's. But soon, they'll publish a better version and I'll use that one for the publication. However, I'll still have a situation of one official version + personal version. As far as I understand, REMAP only works with official versions.

A relatively convenient method would be to create a .chain file, as used by tools like liftOver and crossmap (and probably remap from NCBI). For your use-case, you can simply use the DoSameSpeciesLiftOver.pl command as described in the link. The resulting chain files are pretty self-explanatory and can be used to map BED files of regions of interest.

I am following all the steps. I am stuck here: faToTwoBit ../../genomes/genbank/GCA_000004515.3_Glycine_max_v2.0_genomic.fna.gz GCA_000004515.3_Glycine_max_v2.0.2bit
dyld: Library not loaded: /usr/local/opt/mysql/lib/libmysqlclient.20.dylib
Referenced from: /usr/local/bin/faToTwoBit
Reason: image not found

I suggest you install that with conda instead, then you'll get the mysqlclient dependency (presumably). BTW, I expect you need to gunzip the fasta file.

I have the following error: "ssh: connect to host localhost port 22: Connection refused" for the command "time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=pwd
-ooc=pwd/${target}.ooc -fileServer=localhost -localTmp="/dev/shm"
-bigClusterHub=localhost -dbHost=localhost -workhorse=localhost
-target2Bit=pwd/${target}.2bit -targetSizes=pwd/${target}.chrom.sizes
-query2Bit=../${query}/${query}.2bit
-querySizes=../${query}/${query}.chrom.sizes ${target} ${query}) > do.log 2>&1". I don't understand why it's trying to ssh.

I don’t know where ssh is coming from there either.

This post helped with the ssh error: https://stackoverflow.com/questions/17335728/connect-to-host-localhost-port-22-connection-refused#22255174

Ah, it uses ssh to submit commands to itself as an ersatz cluster.


text:
Where can I find a GWAS data set about genotypes and phenotypes for Zebrafish and rice?
My goal is to work on GWAS method development as a research project. However, human genome data are usually confidential because of the identification problem, so it's very hard to get them. (application takes a long time and I may not get it at the end.)
Therefore, maybe I can try to play with the data of other species first. So I wonder where I can find some GWAS data for Zebrafish and rice. (The reason I am asking for these two is that I find researchers tend to perform genome analysis to mice, Arabidopsis Thaliana, zebrafish, and rice, and I already found where to download the GWAS data for the first two.)
I can find several data sets (like ZFIN) that have the annotated genome information, but that's only one genome. Where can I download the raw data (like genotypes of 500 zebrafishes and corresponding phenotypes) to start my own analysis?

Possible duplicate of Where can I find a database that has phenotype information together with associated SNPs?

@KamilSJaron Thanks. In my defense, the previous question is asking for the associated results that people already find with their own research. This question, I'm looking for a data set that I can research on. What are the differences between human and these species? I did this mainly because it take a long time to apply for human genome data and we may not get it.

You do not need to defend yourself, I am trying to help you in making question clearer. I would recommend you to [edit] your question to add there more details (including those you wrote in comment). Tell us what you want, what you have tried and why you think it did not work.
If the problem with human is dataset size, why not just subsample human dataset? However, GWAS is notoriously known for requirement of huge sample sizes.

I think it will be more productive for you if you choose one of the public human datasets and practise your GWAS methods on that. A couple of examples examples:

1000 genomes -- use "population" as your trait of association. You'll find that you get association across the whole genome, so the analysis method would change from "what bits are associated" to "what are the bits that are most consistently associated"
NCI Genomic Data Commons -- open access (and controlled access) data associated with cancer

I think you need bit more precise phenotype information than population, if you want to test a potential of method. For population you can just look as 2. and 3. axis of PCR and it is going to work (you do not need any formal statistics)...

But that's not going to tell you what markers are associated. I've had success in the past doing a GWAS-like analysis on populations.

You can look at the projection of loci to principal components and therefore you see what loci are actually responsible for the separation.
I do not doubt that you can associate loci with populations. I am trying to say that you do not need any sophisticated method for that because the signal is very strong. I had an impression that the challenge of GWAS studies lays in identification of weak signals, therefore for benchmarking taking "easy" data is not really useful (maybe for sanity tests).


text:
What will be an appropriate mathematical distribution for SNP data?
I found that several papers describe SNPs as a binomial distribution with the probability of "success" equals to minor allele frequency.
However, in my experiments, when I generate SNP array following this distribution, the simulation results behave very differently from the SNP array generated following the random mating procedures. I feel like there is some missing information about the description of a binomial distribution.
I wonder what will be a better mathematical description of the SNP array data? Besides a direct answer, any thoughts, or any suggestion of papers are also highly welcome.
EDIT:
Thanks for these suggestions, after reading relevant materials, I think what I observe is about the Hardy-Weinberg equilibrium affected by evolutionary influences like mutation and migration. I will continue to look and list findings here if I find anything interesting. Please feel free to suggest some mathematical descriptions.
Also thanks for the suggestion of tools. They look very helpful, but I really hope to understand some mathematics behind it.

Are you familiar with linkage disequilibrium and Hardy-Weinberg equilibrium? Looking up review papers on those might give more insights.

Re the edit: given the very low mutation rates in eukaryotes, mutation effect on HWE in SNPs should be negligible - or at least, much lower than the effect of genotyping errors.

Simulating genotypes with realistic correlation structures is indeed not so simple, and there's quite a few papers dedicated entirely to that (e.g. https://bmcgenet.biomedcentral.com/articles/10.1186/s12863-015-0173-4). Also, DEPICT (https://data.broadinstitute.org/mpg/depict/index.html) comes with a number of simulated GWASs to generate the nulls, so that's a simple way to obtain some ready-made fake data.
On the other hand, indepedent SNPs should behave as draws from $Binomial(2, MAF)$. Note that the standard QC procedure of filtering out genotypes that deviate from Hardy-Weinberg equilibrium is just a goodness-of-fit test against $Binomial(2, \hat{MAF})$. If you employ this filter, some LD-pruning, and still get different results, consider posting more details about the discrepancies you see - that would be indeed unexpected.

Thanks. I agree that independent SNPs behave as draws from Binomial. I think the discrepancies I observe is mainly on a population level.


text:
How to aggregate multiple MAF files together?
I have cancer genomic data (tumor/normal whole exome sequencing) from 50
patients that received the same type of treatment, half of whom responded. These come in the form of 50 .maf files, along with a supplemental file that, along with some other fields, has the Response (Responder vs Non-Responder) field. My question is how to aggregate all of this together so that I can perform some statistical tests on the data. I have a clever way (taken from this blog post) of reading in the 50 .maf files and combining them, but I wonder if this is an appropriate approach.
sample_info <- readr::read_tsv(file = "path/to/sample_info/sample-information.tsv")
maf_files <- fs::dir_ls("path/to/mafs/")
patients_data <- maf_files %>%
  purrr::map_dfr(read_tsv, col_types = list(Chromosome = col_character()))

My thought then was to just dplyr::left_join() the sample_info with patients_data like
patients_data_final <- patients_data %>% dplyr::left_join(sample_info, by = c("Tumor_Sample_Barcode", "Matched_Norm_Sample_Barcode"))

For clarity, here are the column names of both dataframes
> colnames(patients_data)
 [1] "Hugo_Symbol"                 "Chromosome"                  "Start_position"
 [4] "End_position"                "Variant_Classification"      "Variant_Type"
 [7] "Reference_Allele"            "Tumor_Seq_Allele1"           "Tumor_Seq_Allele2"
[10] "Tumor_Sample_Barcode"        "Matched_Norm_Sample_Barcode" "Protein_Change"
[13] "t_alt_count"                 "t_ref_count"

> colnames(sample_info)
[1] "Patient_ID"                     "Tumor_Sample_Barcode"           "Matched_Norm_Sample_Barcode"
[4] "Response"                       "Silent_mutations_per_Mb"        "Nonsynonymous_mutations_per_Mb"
[7] "Mutations_per_Mb"

My task is to find genes that "are enriched in patients who either
responded or not." So as a supplemental question, if anyone has suggestions on which statistical test to use (or how to go about deciding), I'd appreciate that as well.
PS: I am aware of the maftools package which probably has an easy solution to this, but unfortunately my machine is old (Late 2011 MacBook Pro) and unable to run it. (Old Mac --> Can't update OS --> Cant update version of R --> Can't install necessary packages to run maftools)

Ask someone with a newer computer for help, or fire up an AWS instance with Ubuntu and do it there. You will face a lot of difficulties in bioinformatics if your machine is >5 years old and you don't have access to any HPC. Old notebooks can work well as terminals but not as workhorses.