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.


text:
Contamination on genome assembly
I had a question for the community.
I have a genome of a new species that has been sequenced via 150pb Illumina paired-end.
To verify the quality of the assembly I used the BUSCO pipeline which allows me to search for 99% of the BUSCO genes (gene not duplicated) present in a taxon in the genome of my species.
Globally I only miss 18% of the BUSCO genes, which means that the assembly is not too bad.
Now when I plot the G+C% and the coverage of all the scaffolds containing BUSCO genes (so a priori non-contaminating) I should find a fairly homogeneous cloud.
Instead I end up with two different coverage profiles (noted red and blue below):

I then thought that the low coverage sequencing profile around 4x might be a eukaryotic contaminant like a acarians, however when I blast the BUSCOs present on these low coverage scaffolds I get hits that are very close to the genome of my species...
Does anyone have any idea what these are?
here is an output of the contigsTable.csv where the two first row are scaffold with low coverage in the figure, and two last are scaffold with high coverage in the figure

FRC output

Coverage distribution

Here is an Icarus summary of the assembly :
Assembly
# contigs (>= 0 bp) 1038631
# contigs (>= 1000 bp) 120132
# contigs (>= 5000 bp) 22245
# contigs (>= 10000 bp) 9986
# contigs (>= 25000 bp) 1322
# contigs (>= 50000 bp) 97
Total length (>= 0 bp) 655530043
Total length (>= 1000 bp) 451497861
Total length (>= 5000 bp) 260894112
Total length (>= 10000 bp) 174183054
Total length (>= 25000 bp) 45006012
Total length (>= 50000 bp) 5956145
# contigs 214542
Largest contig 120620
Total length 519138286
GC (%) 32.94
N50 5062
N75 1624
L50 21982
L75 71030
# N's per 100 kbp 22.28

Here is the BUSCO summary as well :
2544 Complete BUSCOs (C)
2312 Complete and single-copy BUSCOs (S)
232 Complete and duplicated BUSCOs (D)
1060 Fragmented BUSCOs (F)
811 Missing BUSCOs (M)
4415 Total BUSCO groups searched
C:57.7%[S:52.4%,D:5.3%],F:24.0%,M:18.3%,n:4415

EDIT
I ran the fastq sequences with another aligner (MEGAHIT) to see if it was not because of an assembly issue and I got the same profils coverage see here :

This could be organism-specific. We don't have a lot of info so far, so I would check a few more things:
Run something like FRC_align. Check if there's a clear signal between regions flagged as suspicious by it and your coverage graph.
Is it a eukaryote? Plant? Check where mitchondria and chloroplasts are on the plot. They will have different GC/coverage signals than the rest of your assembly, which is fine, expected even.
Worried about contamination? Run kraken2 on your raw data and see what you get there.
Run RepeatMasker and check if repetive regions contribute to the effect you see.
Finally, unless your specimen is low heterozygosity, haploid, or inbred, you probably have some effect from ploidy. There are numerous ways to check that. A self-self dot plot with MUMmer and looking coverage histograms are a good start. The idea is that you should see two peaks (let's assume diploid). One for the split haplotigs, and one for collapsed.
I checked the GC distribution of 45 genomes in your order. It's kind of all over the place, but there are some bimodal ones:

What does GC vs Length look like? Here's the one from the 45 genomes mentioned above:

Thank you for the tips I will try them and let you know. It is a genome from a male Hymenoptera (so haploid genome).

I ran FRC_align but to be honest I do not really understand what I have to look into the FRC_align output (they are many). Is there one that gives you suspicious candidats?

Go by OUTPUT_HEADER_Features.txt. It will have features for abnormal coverage as well as abnormal PE alignments - given your input size. Compression is where they align too close together, stretch is the opposite. I also edited my post with some more info on genome GC densities from these genomes.

Thank you very much for that. I added the outputs I got, I do not have OUTPUT_HEADER_Features.txt... But I'm wondering is the issue does not come from the fact that there are 2 coverage profils rather than GC content ? I also added content of the file contigTable.csv to show you 4 scaffolds, maybe you will se what is going on there ?

The OUTPUT_HEADER part is just a placeholder. In your image, the file is Features.txt. I suggest if there are any more concerns of that nature, we should move this to chat, as comments are discouraged from becoming extended discussions.

Yes sure: here is a room where we can talk  https://chat.stackoverflow.com/rooms/224463/chippycentra

That is indeed puzzling. My best guess would be contamination, even if they carry some of the BUSCOs. I run all de-novo genomes through blobtools, which creates a similar plot to the one you made manually, but also adds taxonomical annotations. It can look like this:

Also, what is the proportion of the genome that is low-coverage? Does it have a significant span? I almost always see a bunch of tiny contigs with low coverage in de novo assemblies, but they usually represent a very small fraction of the genome.
-- edit --
One more thought, are you sure you don't have a tetraploid species? The two coverage peaks kind of correspond to 8 and 32x right? One is 1/4 of the other. We developed some tools to investigate this kind of problems in data.
-- edit 2 --
Just a few more comments. Coverage differences can have various reasons, but all the sequencing biases would cause a skew of the distribution, hence the apparent bimodality of your coverage (more apparent from the two 2d plots) suggests two distinct sources. In general, could be due to various ploidy levels or contamination either by bacteria, endosymbionts or even other species getting into DNA isolation. Now that polyploidy is ruled out, it boils down to contamination, and to figure out what contamination it is, I recommend the tool I liked above.

Hello, I have 553 265 scaffolds with coverage >0 in total in which 283 812 have covdepth between [0.1-15] (so more than half of scaffolds), in fact when I plot the coverage I get the plot I added on the post.

That sounds like a lot. Are the scaffolds long? What is the sum of the 283 812 scaffolds?

I also added one more thought to the answer. Is it possible you deal with a polyploid?

I don't have much experience but my common sense tells me I would focus on getting larger scaffolds with higher coverage before trying to make fireworks with the content of these scaffolds. If you can't do that for all the genome maybe you should focus on a region and analyze in-depth that region with better data. (addressed to the OP)

In my experience, sequencing data must make sense. If you have extra coverage peak that affects half of your genome, it's a rather big deal and I do think OP is on right track to investigate where it comes from.

I added summary about assembly metrics on the post. And I'm sure it is an haploïd genome since it comes from an Hymenoptera male. I tought it would be contamination but as I said, even if the genome is quite fragmented it still found 82% BUSCO genes. And when I blast the BUSCO in both high and low (two pics) coverage scaffolds, I found an Hymenoptera hit...

Just wanted to show my skepticism on the suitability of methods for assembly evaluation based on evaluating the genes obtained or in the GC percent. Especially with non-robust data. I am not sure how these methods work but seems to me that there are many requirements that could be missed leading to misleading data. But that's an interesting thread and I'm interested to learn.  Please go ahead...

I did the assembly with another assembler (MEGAHIT) and it seems that it does not come from here since I see the same coverage patterns (see figure in edit post)...

Ha, having a hymenopteran male is really useful for ploidy doubts :-). Did you use a single male for the library prep? I thinking if you could have contamination of different lab strains? Is the BUSCO evenly distributed between high and low coverage scaffolds?

source: Stack Exchange

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

sort: votes


text:
How to remove low abundance and less prevalent data from my dataset?
I have a relative abundance table in .tsv format where samples are in column and rows contain the features (pathways). Something like this reproducible example:
# Pathway   ERR2601 ERR26014    ERR260145_Abu   ERR260  ERR2601 ERR2601 ERR2607 ERR2600 ERR26
1CMET2-PWY: N10-formyl-tetrahydrofolate biosynthesis    0.00826631  0.00946472  0.00962571  0.00411153  0.00811879  0.00844421  0.00915602  0.00789428  0.0091798
3-HYDROXYPHENYLACETATE-DEGRADATION-PWY: 4-hydroxyphenylacetate degradation  0   0   0   0.000695608 0   0   0   0   0
7ALPHADEHYDROX-PWY: cholate degradation (bacteria, anaerobic)   0   0   0   0   0   0   0   0   0
AEROBACTINSYN-PWY: aerobactin biosynthesis  0   0   0   0   0   0   0   0   0
ALLANTOINDEG-PWY: superpathway of allantoin degradation in yeast    0.0000331927    0   0.0000335389    0   0   0   0   0.0000189167    0.0000908934
ANAEROFRUCAT-PWY: homolactic fermentation   0.00682776  0.00681683  0.00675732  0.00495716  0.00691969  0.0053589   0.00514634  0.00779249  0.00426403
ANAGLYCOLYSIS-PWY: glycolysis III (from glucose)    0.00866693  0.0100437   0.00829714  0.00524566  0.00960262  0.0075093   0.00716347  0.00916205  0.00866289
ARG+POLYAMINE-SYN: superpathway of arginine and polyamine biosynthesis  0.00106395  0.000544088 0.000774645 0.00370237  0.00102368  0.00133848  0.00148691  0.00114062  0.00188739
ARGDEG-PWY: superpathway of L-arginine, putrescine, and 4-aminobutanoate degradation    0   0   0   0.00164625  0.0000278133    0   0   0   0
ARGININE-SYN4-PWY: L-ornithine de novo biosynthesis 0.00166327  0.00383993  0.000105528 0.00230283  0.00404469  0.00115975  0.000460711 0.00101816  0.00367102
ARGORNPROST-PWY: arginine, ornithine and proline interconversion    0.000132331 0.000244069 0.0000796126    0.000479006 0.000237704 0.000290783 0   0.00066458  0.0000581229
ARGSYN-PWY: L-arginine biosynthesis I (via L-ornithine) 0.0124159   0.00955758  0.0118644   0.00533077  0.00774959  0.0098352   0.0127232   0.00996349  0.0108839

Now, I want to keep features (i.e. pathways) which are present with abundance >0.0001 and present in at least 10% of samples. What bash command can I use in order to do that?

Your "reproducible example" is not reproducible. It should be a tab or comma separated plain text content pasted directly in the post and formatted as code (or uploaded on a GitHub gist). A screenshot of an excel spreadsheet is quite far from reproducibility.

Thanks Ram RS for your response. I shared a google drive link. I didn't know it will be shown as an image. However, I have edited the post and entered the table as code. But, still I don't know why it is not showing as a  .csv table. Please tell me whether you can help me. Thanks for your time.

Your data is a little difficult to understand. Can you describe it better please? Please explain what each row and each column represents. I see duplicate column names (ERR2601) as well as inconsistent prefixes (such as the _Abu prefix on one which seems to indicate some sort of difference from the rest of the columns). If you're looking to pick all rows where at least 10% of columns 2-10 are >0.0001, you should be able to do that in R by counting number of columns that are >0.0001.

Column names are just subject IDs and rownames are the pathway names. So, please don't think about them. Each column denotes the relative abundance of respective pathways. If you add the values of a column, it will give you 1.00. Now, here, as example, I have 9 columns, right? Now, 10% of 9 is 1 (apprx.). So, I want to get an output file containing rows with at least one value > 0.0001.  That means, rows with 1 or more values > 0.0001 will be kept and rows which does not have at least one value >0.0001 will be removed. As, there is not much complex calculation, I would prefer to do in bash.

Bash is the wrong tool for this.

I would look at Python pandas or an R dataframe. The best bash solution would require awk but you'd need to know how to program with it.

Here is a bash/awk solution:
awk -F'\t' 'NR==1{print}; NR>1{N=0;for(i=2; i<(NF+1); i++){if($i >= 0.0001) N++}; if(N>0.1*NF) print}' input.tsv

We make sure the input field separator is tab, then we loop on the number of fields (NF) and add to the tally the N of fields above threshold. At the end we simply print the entire line if N is high enough.
EDITED: for the header, we print the first line (NR==1) as is, then apply the expression when Number of Records > 1.

Thanks a lot @Alexlok. It's working. But, it is also removing the header (i.e. the first line of the table i provided). Is there a way to retain it in the output?

Indeed, my tail command was excluding the header, I edited to apply a different awk rule for the header or the rest of the text.

Thanks again @Alexlok. It's working really good!!!

Assuming you have this stored in a matrix called m, where the pathway labels are row names:
IDX = apply(m, 1, function(x) sum(x > 0.0001)/length(x) >= 0.1)
m_subset = m[which(IDX),]

OP isn't using R. OP also has duplicate colnames.

OP needs to use the right tool for the job, it doesn't matter if they prefer bash.

I agree with this. You'd need to be an @terdon standard coder to be doing this routinely and he'd likely prefer to use awk, in this example. Few coders are that standard. The OP could use awk if they wanted a bash solution.

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


text:
Why SeqIO.parse method isn't working?
I am trying to follow the tutorial from from the Biopython website here  and I am right at the beginning (2.4.1), where I am trying a Simple FASTA parsing example which is
from Bio import SeqIO
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

For an unknown reason, I get the following error:
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-19-c64128d7dcfd> in <module>
      1 from Bio import SeqIO
----> 2 for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
      3     print(seq_record.id)
      4     print(repr(seq_record.seq))
      5     print(len(seq_record))

~\Anaconda3\lib\site-packages\Bio\SeqIO\__init__.py in parse(handle, format, alphabet)
    605     iterator_generator = _FormatToIterator.get(format)
    606     if iterator_generator:
--> 607         return iterator_generator(handle)
    608     if format in AlignIO._FormatToIterator:
    609         # Use Bio.AlignIO to read in the alignments

~\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py in __init__(self, source, alphabet, title2ids)
    181             raise ValueError("The alphabet argument is no longer supported")
    182         self.title2ids = title2ids
--> 183         super().__init__(source, mode="t", fmt="Fasta")
    184
    185     def parse(self, handle):

~\Anaconda3\lib\site-packages\Bio\SeqIO\Interfaces.py in __init__(self, source, alphabet, mode, fmt)
     45             raise ValueError("The alphabet argument is no longer supported")
     46         try:
---> 47             self.stream = open(source, "r" + mode)
     48             self.should_close_stream = True
     49         except TypeError:  # not a path, assume we received a stream

FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'

I am using a Jupyter Notebook.

There's something wrong with your copy of https://raw.githubusercontent.com/biopython/biopython/master/Doc/examples/ls_orchid.fasta

Is there any way not to download a .fasta file every time but say scrape it directly from pubmed genome? This string https://www.ncbi.nlm.nih.gov/nuccore/FJ211590.1?report=fasta isn't accepted as a valid argument. Thank you.

Just download the file

Thanks, this means there's no other way?

If you wanted to program in a cache you could, but that’sd be excessive IMO.

Personally I think its late and we're all tired.

"FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'"

The traceback is pointing to line 2, the SeqIO as you point out.
In my opinion, the script should be in the same directory as the fasta file. You know how to identify the path in bash (pwd or echo $PWD), although you must be using Windows (Anaconda has \ rather than /). Alternatively the file names has a typo. You know the
for seq_record in SeqIO.parse("~\path\ls_orchid.fasta", "fasta"): # where path is the dir(s) leading to ls_orchid.fasta, but obviously use / if its Linux.

I would assume you can alternatively dump the .ipynb file in the location where the fasta file is.
As personal choice I prefer a formal IDE such as 'Visual Studio Code' or "PyCharm". In my personal opinion, I am aware some will disagree, Jupyter notebook good is for running code on a remote commercial server for data science work. However, the traceback - which you clearly have activated - looks good. I can see that you're Anaconda3 installation looks good, you must be using Windows (which I've never worked with) and you're running Jupiter from Anaconda3 ... so cool.

Thanks, Michael, these are my first steps with biopython - thought it would work ok with win10 jupyter nb, and yeah checked the path - it was wrong so now everything's fine. Agree 100% that Linux is better.


text:
Error processing run file with bcl2fastq - Unable to open '.../RunInfo.xml' file for reading
I'm attempting to demux a nextseq run which I've successfully done several times in the past on other runs. I uploaded the run from my PC to the server (Ubuntu 18.04.4 LTS) without issue. Then I ran:
```bcl2fastq --runfolder-dir /data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0```

From which I got the error:

'/data/201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml' file for
reading```

However, the RunInfo.xml file is present and readable:
$ ls -l
total 96
drwxrwxrwx 10 user user 90112 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF
-rwxrwxr-x  1 user user 1929 Nov 11 07:37 SampleSheet.csv
$
$ ls -l 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
-rwxrwxr-x 1 user user 28510 Nov 11 08:08 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
$
$ head 201109_NB551126_0128_AH5G5KBGXF/RunInfo.xml
<?xml version="1.0"?>
<RunInfo xmlns:xsd="http://www.w3.org/2001/XMLSchema" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" Version="4">
  <Run Id="201109_NB551126_0128_AH5G5KBGXF" Number="128">
    <Flowcell>H5G5KBGXF</Flowcell>
    <Instrument>NB551126</Instrument>
    <Date>201109</Date>
    <Reads>
      <Read Number="1" NumCycles="75" IsIndexedRead="N" />
      <Read Number="2" NumCycles="6" IsIndexedRead="Y" />
    </Reads>

Has anyone had this issue or know of a fix? Thanks!
Update: contents of run folder
total 2452
-rwxrwxr-x 1 stewart stewart   977 Nov 11 07:59 '2020-11-09__14_05_58_94F12309 GetBestZ Camera 1 Focus Red top surface.snr.csv'
....
-rwxrwxr-x 1 stewart stewart   244 Nov 11 07:59 '2020-11-09__14_08_54_FE05AB9F SafeParabola called by CalculateSignalStd.snr.csv'
-rwxrwxr-x 1 stewart stewart   253 Nov 11 07:59 '2020-11-09__14_08_54_FE0620BC SafeParabola called by FindPeakUsingParabolaFit.snr.csv'
....
drwxrwxrwx 2 stewart stewart  4096 Nov 11 07:59  Config
drwxrwxrwx 3 stewart stewart  4096 Nov 11 07:59  Data
drwxrwxrwx 3 stewart stewart  4096 Nov 11 08:07  Images
drwxrwxrwx 2 stewart stewart 12288 Nov 11 08:07  InstrumentAnalyticsLogs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  InterOp
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Logs
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  Recipe
-rwxrwxr-x 1 stewart stewart    48 Nov 11 08:08  RTAComplete.txt
-rwxrwxr-x 1 stewart stewart  6200 Nov 11 08:08  RTAConfiguration.xml
drwxrwxrwx 2 stewart stewart  4096 Nov 11 08:08  RTALogs
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead1Complete.txt
-rwxrwxr-x 1 stewart stewart    36 Nov 11 08:08  RTARead2Complete.txt
-rwxrwxr-x 1 stewart stewart   925 Nov 11 08:08  RunCompletionStatus.xml
-rwxrwxr-x 1 stewart stewart 28510 Nov 11 08:08  RunInfo.xml
-rwxrwxr-x 1 stewart stewart 26802 Nov 11 08:08  RunParameters.xml
-rwxrwxr-x 1 stewart stewart  1929 Nov 11 08:08  SampleSheet.csv```

What does ls -l /data/201109_NB551126_0128_AH5G5KBGXF/ give?

It was too long for a comment so see update above. Interestingly there are a lot of .csv files I'm not used to seeing but the files that should be there are present as far as I can tell!

I've got it .. just a minute

I repro the error.
ls /ml

ls: cannot access '/ml': No such file or directory

If run
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

I suspect it will work.
At present you are attempting to call a directory /data, which doesn't exit. If you type ls /etc that will work. The difference between ./data and /data in a file structure is alot, but on a command line it appears minimal.

Wait ... if that doesn't work .. because I suspect the calling directory is data then the following will,
bcl2fastq --runfolder-dir ./201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

OR
bcl2fastq --runfolder-dir ./data/201109_NB551126_0128_AH5G5KBGXF --sample-sheet ./data/SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0

It depends on the path where you are making the call. You could sort this out using $PWD or just echo $PWD then copy and paste the path name.

Such a simple fix! Thanks for your solution. I simply ran this from the folder containing the folder and it worked: bcl2fastq --runfolder-dir 201109_NB551126_0128_AH5G5KBGXF --sample-sheet SampleSheet.csv  --no-lane-splitting --barcode-mismatches 0


text:
two aligners combined results
I would like to ask if someone uses two different aligners to produce count matrices and then run the same alogrithm for DEG analysis, would it make sense to find the intersection of the DEGs in order to conclude to the "real" DEGS?

Never heard of that. Main problem in RNA-seq I see is underpowered experiments, therefore unreliable statistics. Just because one aligner may and the other ma not assign certain reads to genes the intersection does not necessarily increase confidence, it also increases the false-negative rate. I would stick to standard procedures, take the aligner (or pseudo/selective aligner, kallisto or salmon) you like most and then proceed with established statistical pacakges.

Thank you so much for your feedback!Yes that particular problem I am trying to solve..!

Any particular reason you try to come up with that? Is the gene you are interested in not showing the behaviour your expect so now you try to tweak results?

No, I am just trying to understand which algorithm to trust and which approach would give the most reliable results.. I am not trying to tweak results..

I suggest you read papers that actually benchmark aligners, there is much work already done on this. Currently, the field seems to prefer the pseudo- and selective aligners such as kallisto and (my personal favourit) salmon. I would always read existing benchmarks first before trying to come up with personal solutions.

This sort of procedure is done in tool papers to benchmark results compared to different aligners/mappers. See the salmon or kallisto papers for examples. Using the intersect of DE results from different mappers will end up yielding somewhat conservative results, not unlike using different differential expression packages and taking the intersect of their outputs. In practice you could just take a more conservative p-value cut-off and save the hassle of using multiple tools, the results would end up being similar.

Yes, I have read in tutorials that the best approach would be to find the intersection of DEGs produced by different algorithms (deseq2, limmavoom, edger) and that s the reason why I wondered if I should the same with the aligner as well.. Thank you for your answer!


text:
Opening gctx files
I installed cmapPy but I don't know how to use it to read gctx files. I read the documentation. If anyone knows how to do it with cmapPy or any faster method your help would be greatly appreciated. I don't know a lot about this stuff so sorry for the nonspecificity of the question.

Hi @Bluebells you question is perfectly fine. Could you explain your project please?

@Michael There is no project, I'm not even an academic otherwise I could just use clue.io. I'm just interested in seeing how results from GWA studies might be used for drug repurposing. I've been using cmap in conjunction with predixcan. I've been experimenting a lot with gene expression signatures and how they might translate into therapeutic candidates. Up until now I was using a previous dataset on broadinstitute.org/cmap but the site is down, probably forever.It's a hobby I guess.

Okay thats clearer thanks. Using metadata matrix manipulation interfaced with GWA. Its written in Python2 ... that is a headache and its not a good sign.


text:
How to find the number of contigs produced? N50 length in base pairs? N90 length in base pairs?
I am having trouble with some underlying questions about my project.
I have ran a Trinity tool to create a trinity.fasta file.
To determine some underlying questions, I used the utility asm_stats to gather some basic statistics.
I am just having some issues of calculating data.
DATA:
Main genome scaffold N/L50: 7464/1.032 KB
Main genome contig N/L50: 7464/1.032 KB
Main genome scaffold N/L90: 14314/106
Main genome contig N/L90: 14314/106
Number of contigs (scaffold length):
25,534 (ALL)
25,534 (100)
21,609 (250)
10,818 (500)
907 (1 KB)
38 (2.5KB)
5 (5 KB)
2 (10KB)
With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?
Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?
Thanks!

With the data attached above, would my N50 and N90 length in BASE PAIRS so excluding the kb, be 7464 and 14314, respectively?

Looks like this is inverse so your N50 is 1032bp and your N90 is 106bp. N90 cannot be longer than N50 by definition. You L50 and L90 are 7464 and 14314

Would the number of contigs produced by this assembly be 7, since there are 7 different contigs produced or would it be 25,534?

25534
Now for the real answer. Do not calculate or report N50 on transcript assemblies. It is more than useless. Transcriptomes are - by nature - dynamic in length. Longer is not better in this regard. Are your short transcripts lower quality assemblies, probably not, but N50 treats them as such. Please read the Trinity wiki for other, better ways to QC your transcript assembly.

Hi @BastianSchiffthaler, thank you for your help! Ok so basically, one is the count and one is the length. Therefore, would my N50 length be 1032bp or is that my count? Then would 7464 be my length?

Also, say if my N50 length is base pairs was 1032 and I wanted to find out "how many contigs are as long or longer than then the length of 1032; would this simply just be 25,534 - 1032 = 24,502? Thanks! @Bastian Schiffthaler

Let's say you order your contigs from longest to shortest. After that, you calculate the sum of all contig lengths, i.e. the total length of your assembly. Now you go from the longest contig to the shortest until the cumulative sum reaches >=50% of the total length. The length of the contig at that position is the N50, the number of contigs you had to include to reach 50% is the L50. N50 = length = 1032bp. L50 = count = 7464. 25534 − 7464 = 18070 contigs are shorter than your N50.

So then that would simply just mean that 7464 contigs are as long or longer than my N50 @BastianSchiffthaler and 14314 for N90? I apologize for the repetitive questions. I am just starting to learn bioinformatics, so this is all a new language to me. I appreciate all your help!

Yes. That's correct

Thank you very much! @BastianSchiffthaler

Hi @BastianSchiffthaler, I was wondering if I could get your advice for one of my commands? Thanks!

Hi @Bastian Schiffthaler, I was wondering if I could get your advice for one of my commands? Thanks!

Please open a separate question.


text:
How to print alternative allele from VCF with htslib c++
I would like to print the alternative allele from a vcf file which i load in with htslib.
The description of some functions:
https://github.com/samtools/htslib/blob/develop/htslib/vcf.h
I have used bcf1_t objects by creating a header and reading the information.
htsFile *test_vcf = NULL;

  // creates header
  bcf_hdr_t *test_header = NULL;

  // initialize and allocate bcf1_t object
  bcf1_t *test_record = bcf_init();

  test_vcf = vcf_open("Testfile.vcf", "r");

  // returning a bcf_hdr_t struct
  test_header = bcf_hdr_read(test_vcf);

  while(bcf_read(test_vcf, test_header, test_record) == 0){
    std::cout << "chr" << bcf_hdr_id2name(test_header, test_record->rid) << ":" <<test_record->pos+1 << std::endl;
  }
  bcf_hdr_destroy(test_header);
  bcf_destroy(test_record);
  bcf_close(test_vcf);

But im not able to extract the reference or alternative allele information using this structure bcf1_t.
As far as i can read its perhaps possible with the structure and function:
bcf_dec_t
bcf_unpack()

But i am not able to get it to work since all the other functions seems limited to the bcf1_t struct.
Can anyone help with extracting the REF and ALT sequence?
Thanks for your time

see the code in https://github.com/samtools/htslib/blob/develop/vcf.c#L3104
 bcf_unpack((bcf1_t*)v, BCF_UN_ALL);// if you need genotypes too, might use another constant to avoid decoding all the structure
 (...)
 kputc('\t', s); // ALT
 if (v->n_allele > 1) {
        for (i = 1; i < v->n_allele; ++i) {
            if (i > 1) kputc(',', s);
            kputs(v->d.allele[i], s);
        }
    }

Thanks for reply, i'll try that!.

Okay, so as far as i understand based on this function, s is a kstring_t* which you can create using mem "kstring_t s = (kstring_t)&test_header->mem;", and then all the information is added to the pointer "s". But how are we then able to print or access the information in "s" ? I can see its a pointer declared with a reference, but then i should be able to print the values with *s. Again thanks for your time.

kstring_t is just the implementation of a dynamic string. you could just use cout << v->d.allele[i] << endl;

okay, that was my initial try, i just looked into it, and the problem is my n_allele for some reason is 0. So ill play around with it some more. Thanks again.


text:
How to replace sequence identifiers in a fasta with OTU IDs from another file?
I'm pretty new to Unix and bioinformatics and having a hard time accomplishing the following. I have one FASTA file with sequences and headers, and one OTU ID mapping file (.txt) with OTU IDs and their sequences. The sequence header in the .fasta file partially matches the sequence listed in the OTU ID file. What I want to do is replace the current headers in the FASTA file with just the OTU IDs from the .txt file. The files look like something like this:
OTU_ID.txt
denovo123 U1019.1010.2020_20384
denovo453 U1472.2992.2020_20443
denovo763 U1928.1898.2020_20343

.FASTA
>U1019.1010.2020_20384 0000-AHDJ-38347_223
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
>U1928.1898.2020_20343 0000-AHDJ-24247_242
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

and my desired output is that the .fasta file would look like this:
>denovo123
GAAATACCGAGAGAGAAATACCGAGAGAGAAAGAAATACCGAGAGAGAAATACCGAGAAATACCGAGAGAGAAATACCGAGAGA
>denovo763
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

Is it possible to do this in the command line?

Yep sure we do it all the time. What coding 'languages' do you know? I assume you are wanting a grep based solution.

Thank you! I only know R at the moment, but am starting to work with bash in the command line (I'm using a Mac).

Hmmm was going to suggest a Perl one liner, but extracting the file name (for me at least) is more of a script than a command line exercise.

You can do it with bash quiet easy. First create a backup of your fasta file with
cp otu.fasta otu_backup.fasta

Then you can run:
while IFS=" " read -r otu fastaId rest;
do
sed -E -i "s/>$fastaId .+/>$otu/" otu.fasta; done < ids.txt;
done < OTU_ID.txt

What does the code?
With IFS=" " you declare a whitespace as field seperator and by reading each line of your OTU_ID.txt, you have the two variables otu and fastaId which hold the otu name and the id of the Fasta Header.
Within the while loop, sed is looking each loop for the corresponding id within a RegEx and substitute it with the corresponding otu name. By using the -i option, you overwrite the existing file. This is why you need to create a backup.

You can put you fasta file in a single line per sequence:
perl -pe 'if(/\>/){s/\n/\t/; s/>//; s/\n//}' file.fasta > file.tsv

This line will create the next output:
U1019.1010.2020_20384 0000-AHDJ-38347_223   GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
U1928.1898.2020_20343 0000-AHDJ-24247_242   GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

Then with awk you can match columns separated with tab:
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1]"\n"$3}' otu_id.txt file.tsv

You will get this output:
>denovo123
GAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA
>denovo763
AAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGAGAAATACCGAGAGA

You can add $1 to check if you are getting the right id.
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1], $1, "\n"$3}' otu_id.txt file.tsv

If everything is alright you can redirect to a new file:
awk -F' ' 'FNR==NR{a[$2]=$1;next} ($1 in a) {print ">"a[$1]"\n"$3}' otu_id.txt file.tsv > file_ids.fasta

How it works, well awk implicitly loops through each file, one line at a time.
FNR==NR{a[$2]=$1;next}
NR is the number of lines that awk has read so far and FNR is the number of lines that awk has read from the current file. Thus, if FNR==NR, we are still reading the first-named file: otu_id.txt. For every line in file.tsv, we assign a[$2]=$1.
Here, a is an associative array, and a[$2]=$1 means saving file.tsv first column, denoted $1, as a value in array a using file.tsv second column, $2, as the key.
next tells awk to skip the rest of the commands and start over with the next line.
($1 in a) {print ">"a[$1]"\n"$3}
This means that we are reading the second file: otu_id.tsv. If match the first field of the line in file.tsv, as determined by the contents of array a, it will print out a line with the values of fields a[$1] and $3 separated with a new line character.


text:
Assessing PyMol sequence alignment object
I use cealign in PyMol for structure alignment. Instead of visualization, I want to return the alignment object to my python script for further analysis. Is there a function to return the object without I/O to a clustalw file?
edit:
By alignment object, I meant the sequence alignment in cealign.

Thats a good question. The answer should be 'yes' outside Pymol based on my own experience in Biopython PDB, but in Pymol I don't know. Someone will definitely know around here.

By alignment object do you mean the 4x4 rototranslation matrix applied to the mobile element or do you mean a MSA object? with PyMOL, using Python you can get both painfully (i.e. using transform=0, object ='rototrans' attributes on cealign). For the former, 3D operations are an utter nightmare (the Mad Hatter's tea party in Alice in Wonderland is a mockery of quaternion rotations... and not unduly!) —is this to move a ED map with gemmi? Because there are less maddenning ways.
For the latter, http://tcoffee.crg.cat/apps/tcoffee/do:expresso is probably a better choice.

@MatteoFerla good remark. Let me edit the question. The alignment object I meant is the sequence alignment. Since I am performing RMSD analysis, MSA in PyMol would be more convenient.

In the application, unfortunately the only option is doing

cealign ,, object =
Menu > Export Alignment... on your cealign  object
and converting the Clustal format using one of the millions of tools online

get_fastastr
As you know cealign adds gaps, but cmd.get_fastastr('') strips them. There is no option to keep the gaps —annoyingly especially with missing loop info!
The help has this to say:
 get_fastastr(selection='all', state=-1, quiet=1, key='', _self=<module 'pymol.cmd' from '/Users/matteo/miniconda3/lib/python3.7/site-packages/pymol/cmd.py'>)

Correction and dissection

So I had a quick check to see if there is an alternative route to get the sequence string with gaps. These are my observations using the Python3 module, sorry if they aren't polished up.

In the comments I said something incorrect. cmd.cealign returns a dictionary with the 4x4 rototranslational matrix —no need for an object. Using Python3 and two TIM barrels as an example:
import pymol2

with pymol2.PyMOL() as pymol:
    pymol.cmd.fetch('1TIM', 'tim')
    pymol.cmd.remove('not chain A')
    pymol.cmd.fetch('5AHE', 'hisA')
    x = pymol.cmd.cealign('tim', 'hisA')
    print(x)
# {'alignment_length': 184, 'RMSD': 3.8800425577361755, 'rotation_matrix': [...]}

While the object does not contain any useful info from PyMol directly*:
    ...
    pymol.cmd.cealign('tim', 'hisA', object ='rototrans')
    pymol.cmd.get_session('rototrans')['names'][0][5]
    # useful info is generally in ['names'][0][5]

(*) The code is actually wrapped C++ code, so the object rototrans is in the C++ layer and can only be access by session or possibly by some functions in cmd.exporting. In fact, it can be exported to an aln string (ClustalW).
    aln = pymol.cmd.exporting.get_alnstr('rototrans')

In the PyMOL application, doing the cealign with an object (cmd.cealign('tim', 'hisA', object ='rototrans'), allows the  Menu>Export Alignment...>'rototrans' which is saved as a ClustalW file, as you mention. Namely:
CLUSTAL

tim          APRKFFVGGN----WK-------MNGKRKSLGELIHTLDGAKLSADTEVVCGAPSI------Y
hisA         MIIPALDLIDGTVVRLHQGDYARQRDYGN---DPLPRLQDYAAQGAGVLHLVDLTGAKDPAKR
                           ..       ......   .....*..................      .

tim          LDFARQK---LDAKIGVAAQNCYKVPKGAFTGEISPAMIKDI-GAAWVILGHSERRHVFGESD
hisA         QIPLIKTLVAGVNVPVQVGGGVRTEE--------DVAALLKAGVARVVIGS---------TAV
             .......   ................        ..*..... .*..**..         ...

tim          ELIGQKVAHALAEG---LGVIACIGEKLDEREAGIT-----------EKVVFQETKAIADNVK
hisA         KSPDV-VKGWFERFGAQALVLALDV---RIDEHGTKQVAVSGWQENSGVSLEQLVETYLPVGL
             ..... *.......   ..*.*...   ...*.*..           .....*..........

tim          DWSKVVLAYEPVWAIGTGKTATP------QQAQEVHEKLRGWLKTHVSDAVAVQSRIIYGGSV
hisA         ----KHVLCTDI----------SRDGTLAGSNVSLYEEVCARYP---------QIAFQSSGGI
                 ........          .      .......*.......         *......*..

tim          TGGNCKELASQHDVDGFLVGGASLKPEFVDIINAKH------------
hisA         GDIDDIAALRGTGVRGVIVGRALLEGKFTVKEAIQ?WQNVKGHHHHHH
             .............*.*..**.*.*...*..

However, there is no method in cmd.exporting that gives out just the sequence with gaps.
Whereas the cmd.exporting has the following... none applicable.
 print(dir(pymol.cmd.exporting))
['DEFAULT_ERROR', 'DEFAULT_SUCCESS', 'QuietException', 'SDF', 'SDFRec', 'Shortcut', '__builtins__', '__cached__', '__doc__', '__file__', '__loader__', '__name__', '__package__', '__spec__', '_cmd', '_feedback', '_get_dump_str', '_get_mtl_obj', '_raising', '_resn_to_aa', '_session_convert_legacy', '_unit2px', 'absolute_import', 'assign_atom_types', 'cPickle', 'cache', 'cache_action_dict', 'cache_action_sc', 'cmd', 'copy', 'copy_image', 'export_coords', 'fb_mask', 'fb_module', 'get_alnstr', 'get_bytes', 'get_ccp4str', 'get_cifstr', 'get_fastastr', 'get_maestr', 'get_mol2str', 'get_pdbstr', 'get_pqrstr', 'get_psestr', 'get_sdfstr', 'get_session', 'get_str', 'get_xyzstr', 'io', 'is_dict', 'is_error', 'is_list', 'is_ok', 'is_tuple', 'loadable', 'lock', 'multifilenamegen', 'multifilesave', 'multisave', 'os', 'png', 'print_function', 'pymol', 're', 'save', 'savefunctions', 'selector', 'sys', 'thread', 'traceback', 'unlock']

what does cmd.get_raw_alignment(string name) returns ?  https://pymolwiki.org/index.php/Get_raw_alignment

I was not aware of that func


text:
How would we test for significant changes in expression when many of the normalised counts are clustered together?
I've been handed RNA-seq data with a lot of associated covariates.
This data has been put through the DESeq2 package and as a result normalised the data.
One of the transcripts of interest still has non-normal distribution for two genotypes. This is because many of the normalised counts are now clustered around a value of 4.

This is making it tricky to perform statistical tests because in order to perform parametric tests (in this case I want to perform a ANCOVA or multiple linear regression to include the covariates) the data must be normally distributed.
I also don't want to use a non-parametric test because I can't include all the covariates.
Should I just perform parametric tests ignoring this violation or do I have to just perform non-parametric on all my data for a fair comparrison?
Or is there another alternative option I am missing?
Thanks in advance

I personally find it hard to follow here. Please clarify: What is "a value of 4"? Is this referring to read counts? What sequencing platform is this? Is this standard RNA-seq or some kind of custom protocol? What do these plots show? Is this Z-score, and if so based on which values? You asked a question on RNA-seq before which indicated that you are new to these kinds of data. It was suggested there to follow standard pipelines rather than coming up with costom solutions. Did you do that?

The values are normalised counts. Platform is Illumina. This data was put through the DESeq2 package.

Your values of 4 look like they might have been vst transformed. This has a tendency to make non-expressed genes have a value around 4. The vst is in log-space, which you can then see in your plots as they aren't linearly distributed

why don't you use the glm model in DESeq2?

TLDR- just do the analysis in DESeq2 (it can handle multiple covariates)
If I've interpreted your question correctly there seem to be a few confusions here.
First off, for a multiple linear regression aka ANCOVA, one of the assumptions are the residuals are normally distributed - however this is almost never an issue in practice. The distribution of the data is not so important, you make the assumption that the data has come from a population which is continuous and normally distributed that's why you've chosen to do an ANCOVA. But, RNAseq produces count data - not continuous and bounded at zero. You need to be looking at using a Poisson GLM, or negative binomial GLM as done by DESeq2 (lookup overdispersion)
The normalisation done by DESeq is to account for different library sizes,it doesn't transform the data to make it conform to a normal distribution. You'd need to do some form of data transformation to do that, but don't it's unnecessary.
I'd suggest reading the DESeq vignettes and following their standard pipeline.
Hope that is of some help.

Will the poisson GLM or Negative GLM account for the covariates as well?

You can use the negative binomial GLM with covariates yes

Your assessment is perfectly correct. The easiest approach is to perform a transformation to normalise your data and then you are good to go. DeSeq2 has a standard transformation (I think it is called deseq2) and you can check your results via a Q-Q plot and you are looking for an approximately linear relationship.
There are a multitude of transformations, but the transformation proceedure has been heavily investigated and standardised in RNAseq

If normalisation by deseq2 has failed (Q-Q) you simply examine an alternative transformation. However, do not transformation a transformation. If you cannot normalise the data then you are stuck with non-parametric stats in my personal opinion.
However, mixed models can work if you trap the underlying distribution and in a regression will result in a very low residual.

Hi Michael, this data has already been normalised by deseq2 so that's why im unsure how to proceed now to analyse it as the qqplots are not normal :(

Commented on above @Connor


text:
How does Allele Number get determined in VCF files?
I am looking at the publicly available gnomAD v3 VCF file.
The header describes the AN field as: "Total number of alleles in samples."
This suggests to me that it should be a constant integer for each site.
However, as the example query below shows (unfortunately they removed the link to download this file, and you need to use cloud services to download the current files -> https://gnomad.broadinstitute.org/downloads#v3-variants):
bcftools view -r chr1:15274 -H gnomad.genomes.r3.0.sites.vcf.bgz | \
awk -F["\t;"] '{print $1,$2,$4,$5,$68}' -

This outputs (chr, pos, ref, alt, AN_afr):
chr1 15274 A G AN_afr=19388
chr1 15274 A T AN_afr=17866

, i.e. the AN_afr (AN for African ancestry) is variable for the same site for two different SNVs. I've found other examples where there are three different SNVs, each with its own AN_afr. This seems to be the rule rather than the exception for multi-allelic sites, regardless of the FILTER. Can someone explain why this is?
I understand how it can vary between sites according to whether or not genotypes are called. What I don't understand is why it can vary at the same site? If it is indeed the case that it can vary, then what is the rule for picking a subset of samples for computing the allele frequency for one particular SNV?

It would help to have links to the data (https://storage.googleapis.com/gnomad-public/release/3.0/vcf/genomes/gnomad.genomes.r3.0.sites.vcf.bgz) and also an example of the output that you are getting. The answer here (https://github.com/broadinstitute/gnomad-browser/issues/524) suggests you are correct, but I can't actually see what you are seeing.

See also answers https://www.biostars.org/p/360183/, https://www.biostars.org/p/329348/,


text:
FASTA and PDB: How to specify chain?
For proteins that have multiple chains (e.g., 1EMS), is there an easy way to specify which chain I want to use for blastp?
I cannot imagine that I am the first person to have this problem, but so far my only solution has been to write my own script that gets the chain I want from the .pdb file and convert that to a .fasta file.

Just to be clear, are you searching for a way to extract fasta from PDB of a protein of more than one aa sequence? Or getting a specified one...

Either getting a fasta for a single chain from PDB or specifying which chain to use when using blastp.

WEB VERSION OF BLASTP
When you run a blastp query using a fasta file with multiple chains you can choose which chain's results to examine using a drop-down box that appears after the search has run.
TERMINAL BLASTP
The results for each chain will be listed separately.
For ease of searching you can specify the output to be a csv of your desired format (See blastp -help and look up the -outfmt option for more information).
Note that for some proteins, such as the 1EMS you chose for an example, the multiple chains are identical. When you get the fasta file for that specific protein it will only have one entry. However if you really need a fasta file formated so each chain has its own sequence you can get the fasta file for all PDB entries and extract the specific protein_chain you want.

The NCBI BLAST can be used with the PDB DB (which NCBI has). The PDB codes are stored as 4 letter codes underscore chain, e.g. 1GFL_B. The catch is segment identifiers —but generally they are the same peptide so shouldn't be an issue.
You can search specifically the PDB DB in NCBI (not the RCSB PDB) by setting the database to PDB.

I know how to search the PDB with blast, but I am looking for how to specify the query chain.

The _B in the query (ie. the thing you are search with) is the chain of the PDB structure., whereas the 4-letter code is the PDB structure as appears in the RCSB PDB or PDBe. Namely, NCBI does not store the complex as a single file

Yes, I know this. I do not think I explained my problem well, but thank you for trying to answer it. I figured out the answer last night and have added it below.

The problem is that Blast is an algorithm implemented as standalone, within the NCBI, within Uniprot and (badly) within RCSB PDB. So for getting help being specific with questions helps or else you get unhelpful answers —like mine— or worse get your question closed due to apparent lack of research.


text:
Rejected publication
A publication that we have been working on for over 6 months on COVID-19 has been rejected by 2 journals. The objective of our research was to identify effective compounds against the COVID-19 virus by using Molecular Docking studies performed with AutoDock Vina. Which Journal should we submit our manuscript to? Following are the Rejection emails I received from them. What do you think is the reason for the rejection and what should I do to resolve the issue to get my manuscript successfully published.
Rejection Notification from International Journal of Molecular Sciences
"We are writing to inform you that we will not be able to process your
submission further. Submissions sent for peer-review are selected on the
basis of discipline, novelty and general significance, in addition to the
usual criteria for publication in scholarly journals. Therefore, our decision
is not necessarily a reflection of the quality of your work."
Rejection Notification from International Journal of Biological Macromolecules
"Thank you for submitting your manuscript to International Journal of Biological Macromolecules.
I have reviewed your manuscript and I must regretfully inform you that your work presents information that is specific to a selective audience in its nature and interest and is more appropriate for a more specialized journal. The focus of your present work, a computational study, is unfortunately outside the scope of the journal, therefore, unsuitable for publication in IJBIOMAC.
For IJBIOMAC publication, the computational analysis should be used only as a supporting tool to non-computational analysis.
We appreciate you submitting your manuscript to International Journal of Biological Macromolecules and thank you for giving us the opportunity to consider your work."
Thank you for your time.

Sorry to hear this, clearly IJBIOMAC is looking for wet-lab work, rather than singly dry-lab analysis.

I cannot comment on your specific manuscript since I obviously cannot read it but in general the question is whether a study does add an interesting novelty.
The problem with pure in silico studies is that there is essentially no confirmation towards the in vivo relevance. Drug companies simulate and also test many thousands of compounds in a high-throughput fashion on a daily basis so it is often questionable to what extend individual pure in silico studies without in vivo confirmation do actually add to the bigger picture.
That might be the reason your work was rejected but as said this is a rather general comment and I can neither judge content nor quality of your work so please do not feel offended and read it with a grain of salt.


text:
Is is possible to predict ncRNAs from sequence and homology alone?
I'm working with a set of homologous genes (let's call it gene A) from several bacterial species.
I know (from previously published research) that in gene B (a close paralogue of gene A), there is a non-coding, antisense RNA. This antisense RNA regulates gene B by forming a duplex with the gene B mRNA.
I have noticed that, sometimes, an antisense RNA would be predicted for gene A in the same location in the IMG/MER database. This is not predicted for all of the available gene A sequences, only some of them. IMG/MER uses INFERNAL to predict RNAs, I think.
My question is - is there any way to assess how real this is without wet lab data? As it is, it could be that there's no actual antisense RNA in gene A, and this is a spurious prediction due to homology. Or it could be that the antisense RNA is present in all of the gene A sequences, and the RNA prediction software simply missed some of them. Is there a computational way of testing which of these options (if any) is more likely? Thank you for your time!

https://pubmed.ncbi.nlm.nih.gov/17631615/

@Chris_Rands I don't think that software would work for my current problem (even if it wasn't abandonware)... It seems to distinguish between 'coding' and 'non-coding', which it doesn't feel like it'd help when the antisense RNA is inside the gene (so, it is coding in one frame by default). Thank you for trying though!

Is there any conserved secondary structure in this RNA? One could align all known instances of this RNA and see if there is any significant covariation according to R-scape. That could serve as a computational signal that the RNA is real (although wet lab experiments would still be needed).
Another question is whether this RNA is already in Rfam? If so, then the job of making an alignment of known instances is easier.


text:
How do I lift GWAS results to hg38?
There are some questions about lifting between reference builds, e.g. this one. But there doesn't appear to be a question about lifting a GWAS results file to a new reference build (except off-site). GWAS results files are not VCFs or plink files, so they aren't covered by the usual lifting tools (such as CrossMap or LiftOver). How do I lift a GWAS results file to hg38?
For instance, the Lee et al 2018 GWAS on education has this file:
user@computer:/home/user/Lee_2018$ head GWAS_EA.to10K.txt
MarkerName  CHR POS A1  A2  EAF Beta    SE  Pval
rs9859556   3   49455986    G   T   0.6905  -0.02901    0.00151 4.61e-82
rs7623659   3   49414791    T   C   0.3095  0.02899 0.00151 6.05e-82
rs11917431  3   49644012    C   T   0.6973  -0.0292 0.00152 7.40e-82
rs1873625   3   49666964    C   A   0.6973  -0.02916    0.00152 8.10e-82
rs11921590  3   49644193    T   C   0.6973  -0.02916    0.00152 1.17e-81
rs2352974   3   49890613    T   C   0.5 -0.02569    0.00141 6.35e-74
rs7029718   9   23358495    G   A   0.5646  -0.0244 0.00143 7.82e-65
rs7868984   9   23357826    C   T   0.4388  0.02437 0.00143 8.21e-65
rs11793831  9   23362311    T   G   0.4354  0.02418 0.00143 9.34e-64

Yes. LiftOver tools, if they work on VCF data (some only take BED as input) work on GWAS summary statistics. You could use
https://pypi.org/project/liftover/, for example.

Yeah, I figure out the simpler solution is to make a small VCF with an empty case, and then lift that to hg38 using liftover, then merge back the column of new positions.

It turned out to be easier than expected. You figure out which dbSNP version they used. dbSNP is the classification system for snps based on the rsid format. Pretty much every GWAS publishes their results with the rsids for the variants, as done above as well (MarkerName). So you download one of the dbSNP files from the official site. In this case, I am using version 142 since this was used by the 1000 genomes dataset.
After this, it is a matter of joining the tables. This can be done easily in R (you need about 15GB free memory), we load the dbSNP database into memory. Since we don't need all the columns, we skip some of them to save memory:
#library
library(tidyverse)

#cols
vcf_cols = cols_only(
  `#CHROM` = col_character(),
  `POS` = col_number(),
  `ID` = col_character(),
  `REF` = col_character(),
  `ALT` = col_character()
  )

#this is from the 1kg
#ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/other_mapping_resources/
dbsnp142 = readr::read_tsv(
  "ALL_20141222.dbSNP142_human_GRCh38.snps.vcf"),
  comment = "##",
  col_types = vcf_cols
  )

GWAS_EA = read_tsv(GWAS_EA.to10K.txt,
                   col_types = cols(chrpos37 = col_character(),
                                    chrpos38 = col_character()))

#make chrpos
dbsnp142 %<>% mutate(
  chrpos38 = str_c((CHROM %>% str_replace("chr", "")), ":", POS)
)

#count overlap
(GWAS_EA$MarkerName %in% dbsnp142$ID) %>% table2()

#mutate and merge
GWAS_EA = GWAS_EA %>%
  mutate(
    chrpos37 = str_c(CHR, ":", POS)
  ) %>%
  left_join(dbsnp142 %>% select(ID, chrpos38), by = c("MarkerName" = "ID"))

#update file on disk
GWAS_EA %>% write_tsv("GWAS_EA.to10K.txt")

The two new columns, chrpos37 and chrpos38, have the chr:pos variant names in both versions.

I would advise against using dbSNP rsIDs to perform a liftover as they get updated over time. I recently wrote a software to liftover summary statistics and the documentation includes instructions that show you how to easily generate educational attainment summary statistics for GRCh38

Your answer could be improved with additional supporting information. Please [edit] to add further details, such as citations or documentation, so that others can confirm that your answer is correct. You can find more information on how to write good answers in the help center.

Once the plugins have been correctly installed, you could lift your GWAS results with a command like this:
bcftools +munge -f human_g1k_v37.fasta -C colheaders.tsv GWAS_EA.to10K.txt | bcftools +liftover -- -s human_g1k_v37.fasta -f Homo_sapiens_assembly38.fasta -c hg19ToHg38.over.chain.gz


text:
Detecting Allelic Imbalance
I am interested in detecting allelic imbalance for a credible set of SNPs. Currently, I am looking for some packages, hopefully, for Python, that can assist with this. I am aware of this R package, https://bioconductor.org/packages/3.13/bioc/html/AllelicImbalance.html, but I am hoping someone knows of a similar Python library.
Cheers

You're comparing availability on Bioconductor to availability on Python. BioC is specializing on biological software, which Python does not (or at least doesn't have an equivalent community).

Why not just use R? Jumping between that and python is the norm.

What do you mean by "allelic imbalance". This maybe a term used in RNAseq because everyone seems to understand it, but this is not  a term used in population genetics.

@RamRS, I am more familiar with Python, and dont have alot of experience with R. With that being said, I should probably familiarize myself with R going forward.

It's not about the language itself, @Krizbomb. Python doesn't have a BioC equivalent, AFAIK.

Hi @RamRS, Yes I know that there is not a python equivalent, thanks anyways.


text:
DNA sequence error annotation
Suppose I generate a DNA sequence of the following pattern.
AAGTC
And after being passed through a channel with insertion, deletion and substitution errors obtain the following sequence
AAAGGC
Where there was an insertion error in the first position, second position and deletion in the third position and no error in the final position.
I currently know the error pattern but is there a way to automatically identify the error patterns and annotate them?

Hi and welcome. Could you please edit your post to be clearer, speciically are you obtaining errors in the allele calling? What data pipelines are you using?

Stated differently, "I want to perform pairwise alignment and get the resulting CIGAR string." The CIGAR string is the error pattern (note that there may be a number of equal possibilities).

Can you tell us anything about the data that you want to use to identify such patterns?

The original message is a sequence of length 150. After being passed through a channel with insertion, deletion and substitution errors(or no error) it could have slightly different length.


text:
Discard Candidates if N >= 40%
I have obtained ( interspersed repeats ) libraries in the Palm plant genome. Still, unfortunately, the candidates obtained seems false positive a lot of tandem repeats found, due to denovo identifications tools I assume, then employing TRF tandem repeat finder to mask any tandem repeat, I have succeeded to mask them.
But my question here, how can I discard sequence with a belonging header, if the Candidates masked more than 40% of the full length for each of the identified candidates?
Thank you

Hi @BioInfo I understand the problem, its a good question. Clearly the problem is species specific (viruses don't have tandem repeats). Could you hint as to the species please? It would be useful to know the assembler?

@Michael question updated, I'm working with palm genome.

Thanks better, not perfect, but –Errhhh wait ... what is a "Freut plant" ... happy plant (German), "fruit plant"? I'm no wiser.


text:
Convert a fasta file to Newick format
My problem is (should be) pretty simple.
I have a series of 10 sequences which I wish to convert into a Newick format file - what would be good software to do this? I've tried using Jalview and am getting nowhere.
Any help appreciated.

I am not sure if I follow. Do you want to create a phylogenetic tree?

The best would be Biopython using AlignIO, given they are all in a single file otherise SeqIO. The alignment approach,
from Bio import AlignIO
align = AlignIO.read("/home/path/myseq.fa", "fasta")
print(align)
with open("/home/path/example.nex", "w") as handle:
    SeqIO.write(align, handle, "nexus")

(not debugged)
MEGAX or ClustalX which are nice GUI's that will also work. If you've not got Biopython installed then ClustalX is probably the easiest.

Thanks for your input - I'm a Python novice and don't have Biopython and looked at ClustalX - and then I found a solution using Clustalw which works quite nicely - see my answer - but I gave you a +1 for your contribution which led me to my answer.

@M__'s answer led me to clustal.org and some searching led me to this page (Christoph's).
I tried compiling Clustalx (GUI - problem with qmake) which didn't work, but Clustalw (compiled to exe clustalw2) - the command line tool - worked a treat.
So, I had my clustalw2 executable and my 25 Drosophila alcohol dehydrogenase protein sequences and I ran the following commands (from Christoph's page):
./clustalw2 -align -type=protein -infile=droso.fasta

This produces a droso.aln file which looks like this (only a few lines shown):
D_paulistorum               ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_insularis                 ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_equinoxialis              ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
D_willistoni                ETRVAEKLLEHPTQTTLACAQNFVKAIELNKNGAIWKLDLGTLEPIEWTK
                            *. *:**** ** *.:  *.:*********:*****:** .***.* *::

Then on Christoph's page, there is a section:

Generate a Newick tree file

So, I followed the instructions and ran:
./clustalw2 -tree -infile=droso.aln

and obtained a file droso.ph which looks like a Newick format file to me (some lines):
(
(
(
D_madeirensis:0.00000,
D_subobscura:0.00787)
:0.00257,
D_guanche:0.00530)
:0.00336,
(
(
D_pseudoobscura:0.00000,
D_pseudoobscura_pseudo:0.00000)
:0.00262,
D_persimilis:0.00131)
:0.01239)
:0.01132,
D_ambigua:0.01821)
:0.00839,

So, it looks like I have a path to doing what I wanted in the first place - getting from .fasta to Newick format.
Then, using the programs jalview and figtree, I was able to draw trees and the like!

I think you can do this in Jalview.
In your alignment window, open the Calculate menu and select Calculate Tree or PCA....
This will open a Choose Calculation window which you can use to select what kind of tree to build.
When you click on the Calculate button you will get another window with the tree in it.
In that window you can open the File menu and then select Save as... → Newick Format.  Hopefully this will do what you need.


text:
What is the meaning of letter 'c' in tblastn results before genomic sequence position?
An example of a result. So far I have been unable to find the answer on my own.
">CP000000.0:c2005000-2006000 Escherichia coli strain XYZ chromosome
ATG..."

Please provide more details. Without context this could mean anything or nothing.

What more context should I give? I have the amino acid sequence of a protein and I want to find its genomic sequence in all available strains of a certain bacterium. I have tblastn hit results in .txt format. As I understand  after the accession number there is a colon and after it numbers signifying genomic position of hit sequence. In some cases, however, after the colon there is the letter 'c' and I have no clue what that means.

The c is from the database itself, so one would need to know exactly WHAT you're blasting against.


text:
What is the best way to drop alignment artifacts in sequence objects in biopython?
I have a bunch of sequences created with Roary which I want to preprocess. When I read the sequences (genomics data) I want to translate them to protein sequences using:
import pandas as pd
from Bio.Seq import Seq
from Bio import SeqIO

def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.translate())} for rec in fasta_sequences])
    return sequences

The problem is that .translate() fails with:

TranslationError: Codon 'T--' is invalid

Is there a build-in function to ignore the dashes?

That is some cool comprehension within a dataframe ... of course you need import pandas as pd

I think the .ungap() is what I was looking for.
def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.ungap().translate())} for rec in fasta_sequences])
    return sequences


text:
When to use Ariba vs Roary or other genomics tools to call drug-resistance variants?
I have a large dataset of assemblies of bacterial S. aureus genomes (~1000, same species). I am looking at different tools to estimate the presence of drug-resistance factors such as blaZ or ampC. Though, using three different approaches, I get three different answers. Some of the tools I use are Ariba or Roary after gene annotation with Prokka.
I wonder what are the most reliable tools used in the field to recognize drug-resistance and virulence factors?

I think you are in as good a position as anyone to know the answer within MRSA. Ariba has the 'genomics clout', but I'm not sure they are you're standard @Soren. They are genomicists and today is data science, whether impact factor currently reflects that or not. Personally I'd use phylogeny as a starting point, but you know more than I, about ML training, parameterisation and testing.


text:
Metagenomics: Identifying most common sequences
I am working on a project and used the following command:
vsearch --derep_fulllength filtered_merged.fa -sizeout -relabel Uniq -output dereplicated_filtered_merged.fa

and got the following output:
87373926 nt in 203453 seqs, min 310, max 480, avg 352
Sorting 100%
10981 unique sequences, avg cluster 2.0, median 1, max 1287
Writing output file 100%

The output had provided me with the data that 10981 unique sequences have been identified. But I cant seem to identify how many reads of the most common sequence were present in the input data.
Any suggestions will be kindly appreciated!

Can you post a sample of the output fasta? According to the VSEARCH docs the abundance information should be contained there (https://github.com/torognes/vsearch/releases/download/v2.15.1/vsearch_manual.pdf).

According to the VSEARCH docs, since you have specified --sizeout your abundances have been written into the FASTA headers:

--sizeout
Take into account the abundance annotations present in the input fasta file (search for the pattern ’[>;]size=integer[;]’ in sequence headers). That option is active by default when rereplicating.
Add abundance annotations to the output fasta file (add the pattern ’;size=integer;’ to sequence headers). If --sizein is specified, each unique sequence receives a new abun- dance value corresponding to its total abundance (sum of the abundances of its occur- rences). If --sizein is not specified, input abundances are set to 1, and each unique sequence receives a new abundance value corresponding to its number of occurrences in the input file.


text:
Question of Padj value (very high, infinite) got from dds
After I got dds from my count matrix by DESeq(), I use results() method to get Padj and log2 fold changes for volcano plotting,
then I found some gene's Padj values are infinite and the P-value is very close to 0;
this makes those dots on my volcano plotting located very top to Y-aix and also confused me:

is this normal or not, such a p-value (very close to 0) and -logPadj is infinite;
how should I treat those data: are they very different expressed between samples and control or just abnormal data which should be taken off?

Here follows the table:

table 1 is dds results, -log10(padj) show inf;
table 2 is count matrix from featureCount
table 3 is size factor of samples
Please help me !

It is not normal to get p-value close to zero, especially if you have not a lot of samples. If I look at your table 2, the mean counts don't tally with the genes in table 1... so can you sort this out

there are some instances, for example you knock out or over-express a gene, then of course you find this insanely differentially expressed..

The SE are really low in this case though, rather than huge LFC. Very consistent gene expression?

For log(0) error add an arbitrarily small number to the 0 p.adj values (e.g.10^-99)

@StupidWolf Thank you StupidWolf and Greg for answer my question. I didn't list all the samples out, and one of the compounds has the opposite effect from the others, like: d has a positive effect on empty control, and a b c has a negative effect on empty control; According to the answers, I think this may cause the so much differentially expressed(but i'm not so sure).

what is being tested in this table? Hey our answers are as good as guesses if you just show a table. It's like reading the oracle from bones

The fold changes for those are huge (on the order of 32x for the largest), so it's unsurprising for the adjusted p-values to be absurdly small.
There's no reason to treat them differently.

Agreed. If you want to get rid of Inf you can add a tiny number to the pvalues before the conversion, e.g. p + .Machine$double.xmin.


text:
Create identicon of a DNA sequence
How can I create an identicon (a visual representation) of a DNA sequence?

Sequenticon is a Python library for generating identicons for DNA sequences. For instance, the sequence ATGGTGCA gets converted into the following icon:

A web interface is also provided at EGF CUBA (Collection of Useful Biological Apps): Render Sequenticons
Disclaimer: I'm the current maintainer of Sequenticon

Hey Peter, like WOW, what a post! and what an answer ... thank you. In case anyone is wondering this has REAL potential. Personally I think there needs to be further development, but like its VERY cool.


text:
How to calculate RNA copies in qPCR using R?
I have these Ct data from qPCR for a series of samples A-H, I want to quantify the RNA copies of my samples using R:
Sample_ID   Replicate1_Ct Replicate2_Ct
A           22.92         24.21
B           24.29         22.16
C           24.29         22.16
D           27.91         23.76
E           21.44         22.81
F           16.50         18.31
G           17.79         16.56
H           17.34         19.09

In addition, I have measured a standard (STD) curve of copy number of the RNA in question, that I want to use as a reference:
STD_Ct      STD_copies      Log10(STD_copies)
33          4.00E+03        3
28          4.00E+04        4
24          4.00E+05        5
20          4.00E+06        6
15          4.00E+07        7
12          4.00E+08        8

It seems
This the original data form:

This is pretty unclear- at a minimum we need some explanation of the data table. e.g. why do rows G and H not have standards when all other rows do? Is each row a different  standard sample? Do you have a loading control reaction for each sample?

The first two columns are the the mean of two replicates, I wanted to make it simple, the standard is only 6 rows, the standard has nothing to do with the samples, it was there just to build to the standard curve to use it for the absolute method @MaximilianPress

I will edit the question to makes this clearer.

I updated the question with the table from excel @MaximilianPress , it should reflect a better idea.

Discussion of latest update moved to a chat. I think there is a numerical error in the way the predictions are being generated, but you seem ok with the results so I'm calling this done.

Now that we've straightened out the data, it seems like the simplest solution is to just linearly interpolate the copy number from your standard curve.
When I read your standard curve (std) and your sample data (ct) into R, I see that the log10(Ct) fits the copy number fairly well, whereas the linear is a bit harder:

So we can fit a model to predict copy number based on Ct on the STD data:
> fit = lm(std$Log10.STD_copies. ~ std$STD_Ct)
> fit

Call:
lm(formula = std$Log10.STD_copies. ~ std$STD_Ct)

Coefficients:
(Intercept)   std$STD_Ct
    10.6847      -0.2357

# we plug those numbers into the linear equation,
# taking the mean of the two replicates for a best estimate.

> rowMeans(ct[,2:3]) * -.2357 + 10.6847
[1] 5.130429 5.210567 5.210567 4.595390 5.469837 6.582341
[7] 6.636553 6.391424

So those are log10(copy number) estimates, of a sort. You can get much more elaborate, of course.

I accepted your answer and upvoted for the efforts, however; I'm not getting the correct results, both intercept (45.29) and slope (-4.2435) are not accurate, I suppose the correct formula (to calculate RNA copies) would be Copy number = 10^((Ct - Intercept)/(Slope))

@user432797 I think that those parameters would come if you fit the model lm(std$STD_Ct ~std$Log10.STD_copies.), which is then intercept 45.26 and slope -4.23. I believe that you simply want to switch your independent and dependent variables to do your prediction.

I checked another software and it shows intercept (45.29) and slope (-4.2435) @MaximilianPress so I will proceed with these numbers.

@user432797 I honestly wouldn't recommend that- I would recommend instead fitting the model as lm(std$Log10.STD_copies. ~ std$STD_Ct ), because you are trying to predict copy number from Ct. If you use those parameters that you mention, those are predicting Ct from copy number, which is the opposite of what you want. Obviously, you are free to do as you like, but that does not seem appropriate to me.

When I apply my data to this model https://www.researchgate.net/profile/Jack_Gallup/post/How_to_calculate_the_copy_no_in_qPCR_of_the_samples_from_their_respective_Ct/attachment/59d625aa79197b8077983fcd/AS%3A318703059570689%401452996138275/download/Copy+number+question+RG+1-16-16+JMG.xls , I get the same results from basic calculation, if you have a reference that I can refer to, it would be great @MaximilianPress

I take it back @MaximilianPress I just checked https://www.xxpresspcr.com/glossary/y-intercept/ and y refers to Ct while x refers to log of copies.

@@MaximilianPress for some reason the number copies with  10.6847      -0.2357 is very high, it doesn't reflect the real expectation, is that ok?

@user432797 it sounds like there might still be some numerical error. Are the values that I've given in my answer implausible? I can't see your results.

I updated the question @MaximilianPress, you see with values that you suggested, I get copies number with negative values and pretty exaggerated (table 2). But when I flip x and y, I get the right and approximate figure (table 1).

Let us continue this discussion in chat.


text:
Creating a tab delimited column
I have a blast file produced. I executed a blast(x) command outputting both "qeseqid" and "sseqid":
QRv313_NP342_d0_h2_l9    YN13213
QRv313_NP9080_d0_h1_l1   YN5345
QRv313_NP123_d0_h1_l7    YN756
QRv313_NP123_d0_h1_l113  YN9768
QRv313_NP654_d0_h2_l6    YN432
QRv313_NP8_d0_h1_l1      YN3242
QRv313_NP756_d0_h1_l2    YN85686

I have written a command in nano within command-line to obtain the following desired output:
NP342    YN13213
NP9080   YN5345
NP123    YN756
NP123    YN9768
NP654    YN432
NP8_d0   YN3242
NP756    YN85686

I have written a nano script to provide me a tab delimited column of my query and subject id. I am just having trouble moving forward from here. I am unsure as to how I would modify my script to provide me with my desired output.
import sys
file_object = open(sys.argv[1])

for my_data in file_object:

  list =  my_data.split("\t")

  print (list [0], list [1])

Is there a way t alter my command so I can receive the desired output?
Any suggestions would be kindly appreciated!

as far as I understand your question, you want to reformat your output file that consists of two columns:
QRv313_NP342_d0_h2_l9    YN13213

to extract the second "field" from your first column, where fields are delimited by _:
NP342    YN13213

In general, your code is close. You could just split again by _:
import sys
file_object = open(sys.argv[1])

for my_data in file_object:
  list =  my_data.strip().split("\t")
  print (list[0].split("_")[1], list[1], sep="\t")

Now when I teach students, I usually tell them to always expect the unexpected, so you might also want to add some minimal error checking, so here's a version that checks some
minimal assumptions:

from sys import argv
from os import path

usage = f"""
usage: {argv[0]} <blast_output>
"""

def main(in_file):
    with open(argv[1], "r") as in_handle:
        for lno, line in enumerate(in_handle):
            line = line.strip()
            cols = line.split("\t")
            if len(cols) != 2:
                raise ValueError(
                    f"Error in line {lno}: Expected two columns,"
                    " found {len(cols)}. Offending line: {cols}"
                )
            query = cols[0].split("_")
            if len(query) != 5:
                raise ValueError(
                    f"Error in line {lno}: Expected five fields"
                    "in query string, found {len(query)}. Offending line {query}"
                )
            print(query[1], cols[1], sep="\t")

if __name__ == "__main__":
    if len(argv) != 2:
        print(usage)
        raise ValueError("Illegal number of arguments")
    if not path.exists(argv[1]):
        raise FileNotFoundError(f"File {argv[1]} could not be found. Check spelling")

    main(argv[1])

And this is what it produces:
>python3 ex.py in.txt
NP342   YN13213
NP9080  YN5345
NP123   YN756
NP123   YN9768
NP654   YN432
NP8     YN3242
NP756   YN85686

Yup, that was a typo. The keyword is meant to be sep="\t". Corrected in the script

Keep track of the ID's you've encountered within python and perform any filtering you want according to that. Basic python programming is outside the scope of this site.


text:
What is a sensible number of gene/observations to explain PCA variance?
I am working with a set of RNASeq dataset. I have about 4000 observations (genes) on 20 samples and plotting a PCA I found the clustering doesn't vary much when I use different number of genes, but the % variance of PC1 and PC2 does. I have ~1300 DEG. But how do you select a sensible cutoff for the number of genes/observations to explain PCA variance? Would it be correct to select the number of genes that are DEG which account for >50% of the variance between PC1 and 2?
Few examples:

All genes: PC1(34%), PC2(14%)

Top1000: PC1(45%), PC2(17%)

Top500: PC1(49%), PC2(19%)

Top50: PC1(55%), PC2(24%)

Top 5: PC1(75%), PC2(24%)

Hi @Ecg ... this is a good question. What we need to understand is why not dump all the data into the PCA - it can certainly handle it. What you appear to be saying is 'how can I re-jig PCA to give a result when I get no result'. If you can edit your question to help here that would be great rather than a priori assume PC1 and PC2 comprise max resolution

Yes, you are right, I am trying to understand where to set up the cut off to dump only relevant data into the PCA, however I do not understand when you say "when I get no result".

By result I mean clear clusters within PCA

S the cluster IS HAPPENING, just the less genes I select, the more clear (more reproducibility between replicates etc), just I am unsure on how many genes I should select.

Hi, yes I'd guessed this. The more data the fuzzier the signal. There's possibly a way foward, its not to reduce the amount of data .. in 'big data' its not cool. I can't remember its implemented though @StupidWolf might know. Hmmj looking at 'all genes' trapping a limited' variance that is weird.

Of course, you can subset however you want. You just have to adjust the message and be clear when communicating such a result. In this case, you have to specify that your PC1 explains 50% of the variance of all DE genes, or the top 1000 most variable genes.You however can no longer claim the PC explains 50% of the dataset.

If you are using DESeq2, you are performing a PCA on the top 1000 most variables genes and plotting that. Some of your DEGs might be among these genes but there might be other factors other than the condition of your interest. Agree with Bastian its a a matter of what you want to convey with this PCA

And it's not very clear from your question..

This is not really what PCA is for, if I understand the question correctly. You don't select observations to "explain" PCA variance, that is more a feature of inferential tools like e.g. ANOVA. I would prefer to use some other measure to discard genes that are simply not informative or measured accurately. For example, first throw out low-expressed or non-expressed genes that are in the noise range. If what you want is to find a minimal set of genes that discriminate between two "clusters", then there are much better tools to do that than PCA, e.g. LDA or something like that.

Thanks for the discussion from all of you. Just to clarify, my PCA shows three main clusters which coincide with my "conditions", thefeore my PCA would be a visual supporting information that the genes are expressed differently between these three clusters. So I believe, like @Devon mentioned below that the number do not really matter, just by showing the top most variable genes should show enough evidence that these clusters are separated, because using all genes both PC1 and 2 already account for a reasonable "big" percentage of variance (48%).

The general goal of PCA in RNA-seq can be stated as, "I'd like a low-dimension representation of my data to allow easy assessment of the gross structure of my samples, specifically for assessing missing batch effects, samples swaps, etc." Ideally, this low-dimension representation can fit into a plot or two. In other words, we'd like the first 2 (maybe 3) PCs to explain a reasonable amount of variation. A few things to note from this alone:

The number of differentially expressed genes were not mentioned. They have no relevance here.
There was no mention made of the treatment groups clustering together. Seeing this is NOT a goal of PCA.

As you have noticed, the more input rows (genes, assuming you haven't already transposed it) the less variation explained by the first 2 principal components (PCs). On the flip side, the more rows you input the better PCA is describing the structure of your entire dataset.
Question: Do I really care about the gross structure of the entire dataset?
Answer: No, I just want to see problems that I might need to account for.
So you really don't have to use all of the data in the PCA, just the bits that will allow you to see any issues. In general that'll be the few hundred most-variable genes (after some reasonable transformation so you're not just looking at the most highly expressed genes!). The exact number shouldn't really matter very much, something in the range of 300-500 most variable genes should suffice. A common practice is to just randomly pick a nice round number somewhere in or around this range. This is not a number that needs optimization unless perhaps you have a very large number of samples and therefore PC1 and PC2 only account for a small (say 10-20%) percentage of the variation, in which case you're better off looking at either a number of difference principal components or using a different dimension reduction technique, such as UMAP.

Thanks, this was a very clear answer. Like you said, when using all the genes, both PC1 and 2 account already for almosth 50% of variance which is "strong" enough to my belief and to support my dataset. Thnaks again for the discussion above.


text:
Why don't the "child" peptide masses sum to the parent in this De novo sequencing example?

In that example, 129+216+327 does not equal the parent's mass of 492.
My understanding is that during MS, the protein is broken down into peptide masses of varying length of approx. 5 to 20 amino acids. They also lose ions such as water or ammonia like in the example., the experimental spectra I thought would show the masses of these "children" that can be summed up to the parent (if you also add all the ammonia and water ions that are broken off). But that is not the case. Where am I wrong in my understanding?

In a peptide tandem mass spectrometry experiment, each fragment (child) ion is typically the result of breaking a single bond in the peptide backbone yield b and y ions. Because each ion fragment ion is the result of one fragmentation event, we don't expect the fragment ion masses to sum to the precursor mass.
If your peptide was fragmented well, you should get a ladder of b and y ions where the differences in masses between ions in the b and y ion series correspond to the masses of amino acids (+/- any modifications and neutral losses, such as the NH3 and H2O mentioned above).
A great resource for understanding these concepts is this tutorial:
https://www.ionsource.com/tutorial/DeNovo/introduction.htm


text:
How do I rename output folder created using WDL's glob function, with cromwell execution engine?
I made a wdl (genomics pipeline language by the Broad Institute) task for fastqc that outputs the zip files, of which there will be as many as reads that the sample I'm working on at that moment has. We create the output using glob:
  output {
    Array[File] fastqcArray = glob("*.zip")
  }

I would like to give the folder a name that is representative of what it contains, for example fastqc, rather than glob-xxxxxxx, so that it's clear what it is in my task's output folder. Is this possible, and if so, how can I achieve it?

Hi @Mar please do edit your post to mention that WDL is a genomics pipeline (Broad Institute). I recognise without knowing this no-one could help but it maximises transparency

Ok, thanks to Ruben Vorderman on the openwdl github page, I found a solution that does not involve renaming the glob folder (which I believe is not possible?). The only way is to write out the names of the files and pass them explicitly to the Array[File] output. This is how I did it (simplified), in case it helps anyone else in the future:
task fastqc_task {
  File forwardReads
  File reverseReads
  String fwdName = sub(basename(forwardReads),".f.*q.*$","")
  String revName = sub(basename(reverseReads),".f.*q.*$","")
  command {
    mkdir -p ./qc/fastqc_zip
    fastqc ${forwardReads} ${reverseReads} --outdir ./qc/fastqc_zip
  }
  output {
    Array[File] fastqcArray = ["./qc/fastqc_zip/${fwdName}_fastqc.zip", "./qc/fastqc_zip/${revName}_fastqc.zip"]
  }
}