text:
How to identify genomic regions / peaks associated with enhancers (TF binding sites)? Is there a tool or a formal recipe?
I am attempting to identify DNA sequences/regions associated with enhancers in my NGS sequencing data. How can I identify (from peak files/counts for chip-seq, atac etc.) genomics regions/peaks associated with enhancers? I know proximity to genes and possibly histone (i.e. h3k27ac) activity near known genes is possibly an indicator, but is this the only way and how is the computational process formalized/captured? What are the specific set of tools involved? i.e. Peak annotations etc.? Thank you.
As I understand it in my limited bio knowledge, enhancers are at varying distances from genes and sometimes they are multiples (super-enhancers). Are all enhancers already "known/described" (akin to genes)? What is/was the process of finding them (human manual procedure)? It's that more manual/heuristic process (without ML) I am looking for so we can consider formalizing and automating it.
Based on a very insightful answer: I think some input/clarity on the following bit would make the answer(s) even more valuable:
Possible list of steps towards this include following tools/functions:

Bowtie (to align fastq to a genome, mouse, human etc.)

Homer/makeTagDirectory (Sep files, associate chr/coords with tags/sequence reads)

Homer/findPeaks (identify and list genomic regions/sequences with count values above a threshold)

Homer/pos2bed (convert peaks to bed file, genome browser compatible form)

Homer/annotatePeaks.pl (Runs on *bed file - Don't know much about this, but intuitively associate a label with peaks found in 4?)

Question: Is there a step 6? Where do we go after annotatedPeaks for example to get to enhancers?

If you have reads from chip-seq or atac you can use tools like homer or macs to identify peaks in your data. Once you have the peaks, these programs allow you to do motif discovery to annotate your peaks with e.g. what transcription factor binding sites they contain.

Thank you @Niek. Very insightful contribution, that has allowed me to add specificity to the question which I have updated with your input.


text:
Does standardised format to write nutrients of a product exist?
Almost every food product has nutritional values labelling. Is there a digital format, one that would be standardised for use in such a case?
For example I want to create a QR code which is going to contain nutritional facts of a product but I want some other person/applications to be able to read this information without having to adjust their parsing code to get the information required. For this I would need to follow some standardised rules to write down this information, so that everyone can easily parse it, thus I would need some rules that describe: "This should be a .json format, it should have this and that written like this.".
Reference to what I mean about nutritional facts:
https://www.fao.org/3/y2770e/y2770e06.htm#TopOfPage
I have tried to search for it and I couldn't find anything. Also I could definitely create my own format if none exists but if there is already one - I bet it's best to use it.
Example:
{
  "nutrients": {
    "Protein": 50,
    "Vitamin A": 800,
    "Vitamin D": 5,
    "Vitamin C": 60,
    "Thiamin": 1.4
  }
}

If you read plain text of this QR code - you will get the same .json text.
Consider this situation:

Similar QR code is printed on milk, fish, meat and other products in a supermarket.
Users can use some kind of application to scan these QR codes.
Each user might use a different application.

In this situation each application would read it and need to parse it to make use of it.
If product companies that package milk, fish, meat, etc. each use a different format to store this information - then each application has to adjust for every possible format the companies provide on their products.
If there would be a standardised format that every company used - then each application would only need to write code to parse it once or even use a predefined parser that everyone shares.
I am looking for such a format and if the standard exists at all.

The reason I am writing this in Bioinformatics is because it's related to computer science and people related to biology are more likely to have encountered this in their experience.

I don't think this is on topic here. More generally the question lacks detail. StackOverFlow might be more suitable.

@M__ Thank you, I have elaborated on the issue.

source: Stack Exchange

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

sort: votes


text:
PCA analysis of samples in a phylogenetic tree
I have a phylogenetic tree.
Each branch ends have samples (s1, s2, ~ s16).

What I want to do is, I want to make PCA analysis plot for each sample.
First I thought each sample has each lineage (expressed in image as (1,1,1,1)~(2,2,2,2)).
So with that coordinates, I can do PCA analysis.
But on the second thought, the difference between (1,1,1,1) and (1,1,1,2) is much less than that of  (1,1,1,1) and (2,1,1,1), but in the analysis, the difference between them will be same.
I am not sure if this kind of approach is suitable for making a PCA plot with a phylogenetic tree.
How can I make this work?

I don't understand. Do you want to somehow represent the phylogenetic information as an ordination plot instead? If that's the case you would be better off going back to your starting data from which the tree was generated, computing distances, and using those for ordination (PCoA probably, or if it's variant data you could just do PCA on that). You don't say anything about having character data for the tips. Maybe this would be clarified by understanding what you mean when you say "sample". Are you trying to replace the phylogeny, or represent something on top of the phylogeny?

What is the tree type, Phylogram or a dendrogram? Your picture suggests a dendrogram. Expressed in image as (1,1,1,1)~(2,2,2,2) ... this is a bit weird, oh I see its a valuation of the nodes. It is not clear what the value '1' represents mutation?

@MaximilianPress The sorting data from which the tree was generated is described here
https://stackoverflow.com/questions/71953348/python-make-newick-format-using-dataframe-with-0s-and-1s
Hope this helps...

@M__ (1,1,1,1) means lineage 1-1-1-1 if you understand this.
I have assigned the samples with the lineage hierarchy represented as 1 or 2.

There appear two basic issues here:

The PCA would represent arbitrary lineage designation not a genetic quantification.
The PCA will weight the first, second ... fourth position equally, i.e. it is not hierarchical. This is very important for the analysis.

Summary The hierarchical -> unweighted PCA will cause analytical problems and should not be used, given it is ultimately arbitrary. If there was a tree specialist reviewing your work it would not be accepted.
The reason is that position one in the hierarchical cluster defines the entire data set - it cuts it in half, whilst position 4 defines a small a very small subtree. PCA will not see the hierarchy it will see position 1 as of equal worth to position 4. Therefore PCA will start associating arbitrary designated position 4 e.g. branch 2, with each other. Thus 1,1,1,2 is a long way from 2,1,1,2 but PCA will miss that it will think the shared '2' means their closer than they really are and find it of equal distance to e.g. 1,2,1,2 - which or course it isn't.
Fixing it You ask how to fix it

You can go back to the original data and use that to create a PCA based on SNPs ... this has been done for other genetic data types, but cannot done reinterpreting a tree. Using this approach there will be a problem with back mutation, but its a parsimony analysis and you would simply state your assumption - ultimately the influence of back mutation concerns the genetic distance between all taxa in the analysis and if thats not big, the error will be small.
Alternatively if you had a phylogram you could construct a genetic matrix and use genetic distance instead. Thats a separate question because its quite complicated.

tSNE is often used to refine a multivariate analysis and I think thats what you getting at, PCA followed by tSNE is valid. A tree isn't any regular multivariate analysis, its a tightly defined approximation within clear mathematical boundaries measuring/reconstructing evolutionary divergence. A tree specialist will see this, whilst a general multivariate stats specialist will miss it.

I didn't expect this accurate and kind answer. I will keep that in mind. Thank you.

Hi @SGKwon you are very welcome, always here to help.

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


text:
Snakemake scatter-gather with wildcard AmbiguousRuleException
My problem is when using Snakemake scatter-gather feature the documentation is basic and i modified my code according to mentioned in this link
rule fastq_fasta:
    input:rules.trimmomatic.output.out_file
    output:"data/trimmed/{sample}.fasta"
    shell:"sed -n '1~4s/^@/>/p;2~4p' {input} > {output}"

rule split:
    input:
        "data/trimmed/{sample}.fasta"
    params:
        scatter_count=config["scatter_count"],
        scatter_item = lambda wildcards: wildcards.scatteritem
    output:
        temp(scatter.split("data/trimmed/{{sample}}_{scatteritem}.fasta"))
    script:
        "scripts/split_files.py"

rule process:
    input:"data/trimmed/{sample}_{scatteritem}.fasta"
    output:"data/processed/{sample}_{scatteritem}.csv"
    script:
        "scripts/process.py"

rule gather:
    input:
        gather.split("data/processed/{{sample}}_{scatteritem}.csv")
    output:
        "data/processed/{sample}.csv"
    shell:
        "cat {input} > {output}"

I added wildcard option but, i got AmbiguousRuleException: Rules fastq_to_fasta(which is previous rule) and split are ambiguous for the file data/trimmed/Ornek_411-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-81-of-8.fasta
I tried lots of things but either rules are not calling or take AmbiguousRuleException. What am i missing, can someone help?

It's not a problem with the scatter-gather feature per-se; Snakemake is just confusing itself with what text should match what wildcards with these filenames.  You could fix this by defining a wildcard constraint for your samples, and/or changing your filename patterns to clarify what should match what.  The docs on scatter-gather show an example with different directory names and file suffixes that prevent this from happening.  Probably the simplest fix would be to follow that idea and use different directory names for the split parts, like trimmed-split and processed-split.
A hint of what's happening is the repeating ...-of-81-... thing, and the mention of ambiguity between rules.  (See also the --debug-dag option.)  If you have a sample named Ornek_411 and you want to make data/processed/Ornek_411.csv, Snakemake will look for data/processed/Ornek_411_1-of-N.csv and so on, and try to figure out how to make those.  But a filename like that could be the output of either the process rule (sample=Ornek_411 and scatteritem=1-of-N) or the gather rule itself again (sample=Ornek_411_1-of-N).  Obviously you want the former and not the latter, but it doesn't know that automatically.

Firstly thanks changing directory name pattern did trick for me. Thanks for explanation of debug-dag option too. i am aware of this logic because i got that before but with just wildcard but i couldn't set scatteritem logic


text:
Multifasta protein alignment
I have around 4000 orthologous proteins for each of 5 species, I'd like to align them and to concatenate alignments. Which tools can I use? And how?
Marco

4000 isn't much. You'll need to start with an aligner. For aligners the common options are in no particular order:

Clustal Omega
MAFFT
FAMSA
Muscle v5

The current one that 'outperforms' the others strictly on a combination of two published benchmarks comprising 'accuracy' and 'speed' is MAFFT. FAMSA isn't relevant here unless you've a really small computer, like an old laptop, it excels at speed.
However, the situation could have changed recently with the  BioRXIV release of Muscle v5. Please note, I have not used it. Here's the Abstract:

On Balifam, ensembles generated by MUSCLE are shown to align an
average of 59% of columns correctly, 13% better than Clustal-omega
(52% correct) and 26% better than MAFFT (47% correct). The ensemble
bootstrap is applied to a previously published tree of RNA viruses ...

On this benchmark they are suggest Clustal Omega is better than MAFFT. Its not been the run of play that is certain MAFFT rose to dominance in recent history. Des Higgins Clustal W and X was historically the aligner of choice by a long shot, then Muscle, then MAFFT ... now its context specific.
Muscle v5 Its difficult to really get to grips with a BioRXIV ms as an algorithm unless you're a specialist in aligner construction. Please note, the 'RNA viruses' is complicated. As an algorithm this author has kudos and has extended that kudos beyond Muscle. The BioRXIV would have been subject to peer-review, but I am not aware Muscle v5 has been accepted for publication.
Historically, Muscle made early advances in accuracy which were duplicated by Clustal Omega at a much later date (2017?).
Essentially, here your choice is between MAFFT and Muscle v5 for this data set (context specific).
concatenate this is less clear but it appears you are performing 4000 alignments with 5 taxa in each? You'll need a data pipeline.
Trimming and concatenation are easily performed in BioPython using BioAlignIO but you would need coding because that can't be done manually. Computationally, BioAlignIO objects are annoyingly expensive - they are SeqIO objectives inside a MultipleAlignment object. However in your case it matters not a jot and would the recommended route all things considered (unless your coding expertise is Java).


text:
Is there a way to extract spike sequence from a Sars-CoV-2 complete genome (preferably in Python)?
I have a complete genome sequence of a Sars-CoV-2 variant, but I am only interested in the Spike sequence?
I should note that I am not a biologist (nor bioinformatics).
The complete genome looks like the following:
AAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAACGAACTTTAAAATCTGTGTGGCTG
TCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAACTAATTACTGTCGTTGACAGGACAC
GAGTAACTCGTCT...

Its length is 30230.
My question is, is there a way to know where exactly the spike sequence lies among this complete genome sequence, and if possible, are there already existing tools to extract it programmatically (in Python).
Thanks in advance.

Spike Wuhan is here
Using Python then use the BioPython Bio.pairwise2 program here
Yours sequence is the reference, the shorter Wuhan query is the query.
The genomic positions position1 and position2 will then be defined via the output of pairwise2, then cut it by using the genomic position -1 as follows:
spike = completegenome[int(position1)-1:int(position2)-1]

-1 is needed because lists start at 0 rather than 1.


text:
Extract single frame with MDAnalysis
How can I do to extract a single frame as a pdb file from a trajectory traj.xtc file, by using MDAnalysis?
I've got a universe
u=mda.Universe('confout.gro','traj.xtc')
Then I selected a residue by doing select_atoms, and I would like to extract a single frame as a pdb file

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

Q to more involved folk in admin: why was Community bot was invoked for this, as it was a case of nobody bothering to answer because it is a boring 'please link me to command X in the documentation' request. I.e. https://docs.mdanalysis.org/stable/documentation_pages/coordinates/init.html#writers

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


text:
PWM (4xM matrices) to Sequence Logo Visualizations in MATLAB or Python
Is there a MATLAB, or else Python (not R) tool for visualizing sequence Logos from probability weight matrices (PWMs)?
This seems like a basic and simple tool to implement, yet just above the threshold in the space of tedious to write my own from scratch - I was surprised I couldn't find a basic function like a
visualizePWM(pwm_matrix)

that would take a pwm_matrix that is a (4 by m) matrix and spit out the visualization of an m-length sequence motif/logo like:

I see all kinds of over-complicated bloated tools that I thought I ask in case I am missing something.

answer from @M__, converted from comments:
This has been answered on stackoverflow using matplotlib for seq logo plot:
https://stackoverflow.com/questions/42615527/sequence-logos-in-matplotlib-aligning-xticks
IMO its only useful for presentation, pretty pictures, thus there's no need for someone to write a library (which is what you're asking). If I'm using a mutation matrix (which I do), I want the raw numbers and significance values I don't want pretty graphs getting in the way.
Conceptually you would move from 2D to 1D, without subplots (which is doable), i.e. its a matrix. My view on code availability is not universal. I've had requests for separate CLIs (useful for someone clearly) ... which just never happened.

The solution proposed in the link is focused primarily on a tool that is captured as the pyseqlogo packet at the end. I have spent a good part of the day trying install the packet and get it to work, but to no avail it's failed with too many errors and complications to list :(.

After much searching and trial and error, located possibly the optimal tool for my question:
seqlogo
https://github.com/betteridiot/seqlogo
Surprisingly not only seqlogo is an elegant and minimal solution in Python (and R) I have discovered it has a relative implementation in MATLAB! It installed and worked the first try - beautiful output with all the output options I needed - rock solid.
Python one liner (not counting the example pwm input generation)
random_ppm = np.random.dirichlet(np.ones(4), size=16)
ppm = seqlogo.Ppm(random_ppm)
seqlogo.seqlogo(ppm  ,format='png',filename='out.png')
Produces the following ('out.png')

And here is MATLAB Port:
https://www.mathworks.com/help/bioinfo/ref/seqlogo.html


text:
how to retain reads with low mapping quality (MAPQ) scores when using samtools view -q
I have been using the -q option of samtools view to filter out reads whose mapping quality (MAPQ) scores are below a given threshold when mapping reads to a reference assembly with either bwa mem or minimap2. However, using samtools idxstats to count total mapped reads and unmapped reads indicates that these reads with lower MAPQ scores are actually being discarded altogether, since the number of unmapped reads reported (last row) is always zero when the BAM file is created with a -q argument, even though fewer reads are being mapped to the assembly.
For example:
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b | samtools sort --reference assembly.fasta > nanopore_mapped_no_q_thresh.bam
samtools idxstats nanopore_mapped_no_q_thresh.bam

gives
contig_1      4751250 9897    0
contig_2      37697   528     0
contig_3      9131    2366    0
contig_4      98997   701     0
contig_5      4539    365     0
contig_6      3933    3164    0
*       0       0       142

while
minimap2 -ax map-ont assembly.fasta reads.fastq | samtools view -b -q 2 | samtools sort --reference assembly.fasta > nanopore_mapped_q2.bam
samtools idxstats nanopore_mapped_q2.bam

gives
contig_1      4751250 9716    0
contig_2      37697   528     0
contig_3      9131    1415    0
contig_4      98997   700     0
contig_5      4539    365     0
contig_6      3933    3161    0
*       0       0       0

Is there an alternative way to filter mapped reads by their MAPQ scores so that the filtered reads are still retained in the BAM file?

Looking at the samtools view docs, I believe that you need to also set the -p flag if you want to retain "filtered" reads in the file and simply have them be unmapped:

The -p, option sets the UNMAP flag on filtered alignments then writes them to the output file.

Yes, this does the trick. This flag was not recognized by the older version of samtools I was using (1.11) but after upgrading to 1.15 it works. Thank you Maximilian!


text:
Is there a way antiSMASH can be used to screen a large number of genomes just for one type of biosynthetic gene cluster?
When we use antismash command with the genome sequence file as follows:
antismash GCA_002095535.1_ASM209553v1_genomic.gbff

it screens the genome for all the possible biosynthetic gene clusters (BGCs).
Is there a way the tool can be directed to screen for only one specific type of BGC (for example, terpene)?

Hi @snap_v I think the answer to the question is about parsing the output. If you present the raw output for a typical search have a go at a bit of code and wrap it up as a new question - I think you'll get a solution.


text:
Reading lines from a fasta format file
file_name = "consen1.txt"
    with open(file_name) as infile:
        new_dict = {}
        for line in infile:
            if line.startswith(">"):
                key = line.strip('\n')
            else:
                line = ''.join(x.strip('\n') for x in infile)
                new_dict[key] = line
    print(new_dict)
    for key,values in new_dict.items():
        gc = values.count('G') + values.count('C')
        content = (gc/len(values)) *100
        new_dict[key] = content
    print(new_dict)

I am reading from a text file, and everything works fine but for some reason it is skipping the line of sequence from the text file. Here's how the text file looks like,
>Rosalind_4478
GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGAC
CGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGC
TTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGAT
CGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCAC
TCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTC
CTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAA
TGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGC
TGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGT
ATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCC
GTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCA
AGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCT
CGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTA
CTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCT
TTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTG
TGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAG
TA

Expected output:
 {'>Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}

Real Output:
{'>Rosalind_4478': 'CGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}

I do realize that I can do this kind of computation easily with modules like biopython, But I wanted to try doing this from scratch.
Thanks in advance for the help.

You can also unlock hundreds of solutions to this problem on the Rosalind website, once you submit the correct answer.

Your approach is a bit risky code ... Anyway not tested but I think this will work:
with open(file_name) as infile:
    new_dict = {}
    for l in infile:
        line = l.rstrip('\n')
        if line.startswith(">"):
            k = line.lstrip('>')
            new_dict[k] = [] # must initiate
            # new_dict = {k:[]} is possible alternative initiation
        else:
            new_dict[key].append(line)

I not sure without testing which dictionary initiation will work. People do use defaultdict
from collections import defaultdict
new = defaultdict(def_value)

This circumvents the initiation issue to a large extent (if the dictionary key ain't correctly initiated it will throw a value error, or some sort of bug). However it isn't actually needed and I find it unnecessary.
What I am certain about is the following should be avoided because your code depends on the exact line by line order between the fasta header and the code you are associating it with. You likely break that association with the code below,
line = ''.join(x.strip('\n') for x in infile) # you are calling the entire infile not the line

BTW
from collections import Count ()
new = Count()

Is a very good approach for the next part of your script, i.e. counting %GC

Thank you so much... Now I see my mistakes... BTW You seem really experienced with bioinfo, if it's not too much to ask... what does one have to do in order to become really good at bioinfo?

Thanks, the following is just my opinion for bioinformatics. Enjoy reading technical computing books, e.g. O'Reilly or Packet (they are exciting, technology, advances, applications), a good IDE with breakpoints certainly helps to start (I use Spyder) and for me what's really important is to visualise an algorithm without seeing the code.

TRY : https://rosalind.info/

needed to edit accepted code to make it works:
file_name = "consen1.txt"

with open(file_name) as infile:
    new_dict = {}
    for l in infile:
        line = l.rstrip('\n')
        if line.startswith(">"):
            k = line.lstrip('>')
            new_dict[k] = str() # must initiate like string !!!!!!!!!!!!!
            # new_dict = {k:[]} #is possible alternative initiation
        else:
            new_dict[k] += line

print(new_dict)

for key,values in new_dict.items():
    # print(key, values)
    gc = values.count('G') + values.count('C')
    content = (gc/len(values)) *100
    new_dict[key] = content

print(new_dict)

output :
{'Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}
{'Rosalind_4478': 49.55654101995565}

works using this snippet too. but feel like is not orthodox python:
file_name = "consen1.txt"

with open(file_name) as infile:
        new_dict = {}
        for line in infile:
            if line.startswith(">"):
                key = line.strip('\n')

                print(key)

                line_2 = None

            else:

                if not line_2:

                    line_2 = ''.join(line.strip('\n'))

                else:

                    line_2 += ''.join(line.strip('\n'))

                # print(line_2)

                new_dict[key] = line_2

print(new_dict)

for key,values in new_dict.items():
    gc = values.count('G') + values.count('C')
    content = (gc/len(values)) *100
    new_dict[key] = content

print(new_dict)

output:
>Rosalind_4478
{'>Rosalind_4478': 'GCTACATCAACGCCGTGCATGGCGTGTAAAGACCCGCATTTTACAGGCAGTGTATACGACCGGTAACACGATTTAGCATTTACTGAGCCCCTGCGTTATTGGAATATGAGGGAACGTTGCTTATGTGCAAGACATTCCACCTCATCCATAGATACTCTCGCGGCAGATTCTGTGACGGATCGGATCGCTTCAAAGCACGGTACGGTATTAGACGAATGGGTTAAAACTCTCTACAAGCACTCGTACAACTTATGCACACCAACTTCAGCGTATAACATGAGAGCCCGGAAAGGCATATTCCTTATCTTAAGGTAGGGAATACGAAAACCCGCGTCGTCCGCCCACCTCGACGGTACGGAATGTTTGAGGGGCATCCGGGTCTGGCAGGAACTGTGTGGACTGGTGCGTCATCAGCTGGGCTGATTGCGCGAGCCAGAATCGCATGCGCTTCGAATTTTAGGAGAGCGTGCCGACTCACGTATATCATTGACATGGTGAGGGTAGCAGAATCTACATTGCAGCCTATGACGCTACACGGCCGTATAATCAGTGTTAGTGCAGTATTGTGGCATTTCCGTTTCGCATGGTTAACTGAGGCCAAGGTTTGTGCCGGGGGACCGGATGGCTACCTGGAAACTAGGATCTGGCCGTAGTGGTTCTCGCGAACCCTAGTTTCAAGCAAACGCCAGCCGGAGACCGAAGTTAAACGCATCACTACTACTCACGAGTCTTGGACCAAAACAGAAAGTCGGATGCACCATGTACGAGCTGACGTTGTCTTTATGTCCAGTTGGATCTGGACAGCACATCCGAGATTTTGAATTAAAACTTCTTTTTGTGTGATGCTTGGAACACCCGCCATACAACAAATAATAGCGCCTTAAACAGCCTACGGCACAGTA'}
{'>Rosalind_4478': 49.55654101995565}

be sure to check if they are both working with your input

Thanks @pippo1980. I like the
new_dict[k] = str()
new_dict[k] += line
approach rather than append. The append approach will work. Making += work is cool without using from collections import defaultdict.
The GC% using count() method directly rather than via Count()


text:
How to adjust by multiple variables using ComBat-Seq?
I am trying to adjust my RNA data using ComBat-Seq  (from sva R package) since I realised that there are 3 batches that I need to adjust for:

Place (2 levels: place 1, place 2)
Library Preparation Date (16 levels - different dates)
Type of tube (2 levels: A, B)

I have 960 samples and around 62000 genes.
In my biological matrix, I have: Age, Sex, Group (cases,controls..) and WBC counts.
biol_mat = model.matrix(~Age + as.factor(Sex) + as.factor(Group) + LYMPH + MONO + NEUT, data=phenotype)

In the tutorial of Combat-Seq appears how to adjust by 1 variable but it doesn't tell you how to adjust by more than 1.
I have seen a lot of posts using combat that the only way is to adjust by 1 variable and then, with those results, adjust again by the 2nd variable and so on.
That would be:
Adjust by library prep.
raw.cts_adjustedLibPrep <- ComBat_seq(counts = raw_cts_matrix, batch=batch_libraryprep, group=NULL, covar_mod = biol_mat)

Adjust by library prep + type of tube.
raw.cts_adjusted_LibPrep_TypeTube <- ComBat_seq(counts = raw.cts_adjustedLibPrep, batch=batch_type_tube, group=NULL, covar_mod = biol_mat)

Adjust by library prep + type of tube + place
 raw.cts_adjusted_LibPrep_TypeTube_Place <- ComBat_seq(counts = raw.cts_adjusted_LibPrep_TypeTube, batch=batch_place, group=NULL, covar_mod = biol_mat)

For the first adjustment (library prep) it takes around 15min, but for the second... it has been running for more than 2 days.. I stopped it and launch it again, changing the adjustment but I am not sure if it will work..
Does anybody have an idea about how to fix the problem?
Thanks in advance
Question already asked here


text:
nested samples in design DEseq2
I'm running a DEG analysis with Deseq2 by specifing the following design that includes 12 samples from 6 pts in 4 different conditions:
coldata <- data.frame("sample"=c(
  "T58", "T86", "T87",
  "D58", "D86", "D87",
  "C62", "C80", "C84",
  "C62", "C80", "C84"
))
coldata$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87",
  "c62", "c80", "c84",
  "c62", "c80", "c84"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal",
  "CLEFT", "CLEFT", "CLEFT",
  "CRIGHT", "CRIGHT", "CRIGHT"
)

I would model the differences among the conditions by considering that the same individuals are nested in different conditions
I tried the following code
ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix),
  colData = coldata,
  design = ~ class+Pts
)

but this produces the famous error

"the model matrix is not full rank..."

If I reduce the previous dataset to
coldata2 <- as.data.frame(c(
  "T58", "T86", "T87",
  "D58", "D86", "D87"
))
coldata2$Pts <- c(
  "m58", "m86", "m87",
  "m58", "m86", "m87"
)
coldata2$class <- c(
  "Treated", "Treated", "Treated",
  "Distal", "Distal", "Distal"
)

then everything goes well and I'm able to create a DESeq object with the following design:
ddsMat <- DESeqDataSetFromMatrix(
  countData = as.matrix(matrix[,1:6]),
  colData = coldata2,
  design = ~ class+Pts
)

Please can you provide to me some suggestion on how to solve the problem?
Thank you

It is actually simple. The first two conditions are nested with the Pt's starting with m and the other two conditions are nested with those PT's starting with c. That means you cannot run this design. The only solution is to subset the data into two chunks and analyze separately.

The DESeq2 documentation explains this error as follows:

While most experimental designs run easily using design formula, some design formulas can cause problems and result in the DESeq function returning an error with the text: “the model matrix is not full rank, so the model cannot be fit as specified.” There are two main reasons for this problem: either one or more columns in the model matrix are linear combinations of other columns, or there are levels of factors or combinations of levels of multiple factors which are missing samples.

In your case, I think the error is popping up because there are too many defined values in the "Pts" variable. DESeq2 doesn't require you to separately define replicates as a variable; it assumes that different data columns are taken from different biological replicates.
Getting the models right is tricky, and needs a deeper discussion about what the goal of analysis should be. My guess is that the following should work, but only if you're not interested in comparing expression results from different patients:
coldata$Pts <- c(
  "m", "m", "m",
  "m", "m", "m",
  "c", "c", "c",
  "c", "c", "c"
)
coldata$class <- c(
  "Treated", "Treated", "Treated",
  "Distal",  "Distal",  "Distal",
  "CLEFT",   "CLEFT",   "CLEFT",
  "CRIGHT",  "CRIGHT",  "CRIGHT"
)

This may not be what you want to do; other ways to deal with this are explained in the DESeq2 documentation on the Bioconductor website:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#model-matrix-not-full-rank

check if Pts and class are of factor class or not by using
is.factor(coldata$Pts) and is.factor(coldata$class)
if it is not then convert it into factor by
coldata$Pts = factor(coldata$Pts)
coldata$class = factor(coldata$class)
Then try to design your matrix (A design matrix is a matrix of values of the grouping variable).
I think this would work as it works with
design <- model.matrix(~0+coldata$Pts+coldata$class) in limma package.
You can also check in this RNA-seq analysis workflow of DESeq2 where they also used nested samples for their design matrix.

Your original was,
ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), # colData = coldata, design = ~ class+Pts)

What happens when you use
ddsMat <- DESeqDataSetFromMatrix(countData = as.matrix(matrix), colData = coldata, design = ~ class+Pts)

Thus # is removed. The suspicion is # is causing the code to fall over. #  is used for comments in R so the presence of this within a line of code will cause the following code colData = coldata, design = ~ class+Pts to not be interpreted because it can't be seen.

THe #was only added by mistake when I wrote the question. In the original code the colData is not commented


text:
My model doesn't work in test dataset

I got the data of stomach cancer from TCGA, I found some different
expression genes between cancer and normal sample.
I put aside my normal sample, and divide my cancer sample into
train(80%) and test(20%) groups.
I did Unicox (with survival package), lasso (with glmnet package)
and multicox ( with survival package) on train group, to find out
which genes are more related to survival.
I use coefficients of each genes (that I got from multicox), to
built my model and calculate the risk score. Like that: (expression
of gene1×coefficient gene1 + expression of gene2×coefficient gene2
+...).
I divided my sample in to high risk and low risk groups, to compare
survival between them, using Kaplan-Meier.

In training group, Kaplan-Meier shows significant difference between high and low risk group. However when I use my model on my test group, it is not significant. Why this problem occurs? I assumed my model is overfit of my train group, but how can I fix it? Why does this happen even when I use lasso?

Do you have sufficient power? Are your hazard rates/ratios very different between train and test? It would be possible to have a model that performs exactly as well as it did in the training data, but that shows up with no significant difference simply due to having a smaller sample size in the test data.

Hi, thank u for your answer, yes at first I thought it is because of small sample size in test group, so I tried diving samples, 70-30, 60-40, and even 50-50. In all of the them pvalue in Kaplan meier was <0.001 in train group, and in non of the it was significant in test group

Hi @FaridehJafari could you kindly format your question to make it clear?

There is no overtraining here, it isn't deep learning and there appears no parameterisation. If parameterisation was occurring on the training set that can and does cause over-tightening.
Essentially, you haven't included a data split for parameterisation, i.e. there is no parameterisation. The theory of a parameterisation is it must be distinct from training (and obviously) testing. It allows parameter tuned without overtightening. Lasso regression will require parameterisation.
I don't know about the survival regressions for machine learning, but lasso/ridge regression are standard approaches.

So just to explain, what you do is have a split for parameters. Thus instead of 80:20 split it could be 70:20:10 split where the 10% is used for parameter optimisation. Bad code risks leakage from the parameter split and that can result in over-tightening.


text:
Simulation software recommendation to confirm pipeline validity
I'm trying to detect convergent evolution in different samples of the same species.
I'm just trying to see if past 20 generations the nucleotide frequency of snps under pressure converges somewhat for all samples.
What i want now is to validate my pipeline using a simulation.
I would like to simulate the evolution of populations and see if the pipeline i built would detect the snps from the reads.
I'm using pool-seq and illumina with reads of around 100base pairs.
Can you recommend me some software and maybe point me to the tutorial?


text:
How to calculate distance matrix from fasta?
I’d like to perform the following, but I am having difficulty.
I want to convert amino acid sequences in multifasta files into matrices (or vectors) with Euclidean distance and Pearson correlation coefficient.
If you know a solution, please let me know. Thanking you in advance.
File1
>BDP27474.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL
>BDP27485.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL

File2
>BDD77411.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFLVLLPL
>BDD77515.1 |surface glycoprotein [Severe acute respiratory syndrome coronavirus 2]
MFVFFVLLPL

Code
library(tidyverse)
library(ape)
library(Biostrings)

naseq1 <- readAAStringSet("File1.fasta")
naseq2 <- readAAStringSet("File2.fasta")

dist.aa(naseq1)
dist.aa(naseq2)

##Error in numeric(n * (n - 1)/2) : invalid 'length' argument

Add info.
iOS 10.15, Rstudio cloud

There's no matrix for file1 because MFVFLVLLPL and MFVFLVLLPL are identical

You are probably looking for Biostrings::stringDist to compare both strings


text:
Most reads no_feature with htseq-count
I have mapped reads to the yeast (R64-1-1) genome with STAR like this:
STAR --runThreadN 44 --runMode alignReads --genomeDir /home/user/Documents/scripts/pyseqtools/index/STAR/R64-1-1 --readFilesIn /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R1_001_val_1.fq.gz /home/user/Documents/analyses/TT-Seq/user/trim/WT_3_R2_001_val_2.fq.gz --readFilesCommand zcat --quantMode TranscriptomeSAM GeneCounts --twopassMode Basic --outSAMunmapped None --outSAMattrRGline ID:WT_3 PU:WT_3 SM:WT_3 LB:unknown PL:illumina --outSAMtype BAM Unsorted --outTmpDir /home/user/Documents/analyses/TT-Seq/user/tmp --outFileNamePrefix /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3

I wanted to count the reads using htseq-count as follows:
htseq-count --quiet --format bam --stranded yes /home/user/Documents/analyses/TT-Seq/user/bam/R64-1-1/WT_3/WT_3Aligned.out.bam /home/user/Documents/scripts/pyseqtools/gtf/R64-1-1/Saccharomyces_cerevisiae.R64-1-1.106.gtf > /home/user/Documents/analyses/TT-Seq/user/htseq-count/WT_3.count.txt

The problem is that most reads map to no_feature:
__no_feature    4462786
__ambiguous 12644
__too_low_aQual 0
__not_aligned   4907
__alignment_not_unique  12529444

I have check whether the chromosome names in the BAM file and the GTF file matched (which they seem to do):
BAM:
@HD VN:1.4
@SQ SN:I    LN:230218
@SQ SN:II   LN:813184
@SQ SN:III  LN:316620
@SQ SN:IV   LN:1531933
@SQ SN:V    LN:576874
@SQ SN:VI   LN:270161
@SQ SN:VII  LN:1090940
@SQ SN:VIII LN:562643
@SQ SN:IX   LN:439888
@SQ SN:X    LN:745751
@SQ SN:XI   LN:666816
@SQ SN:XII  LN:1078177
@SQ SN:XIII LN:924431
@SQ SN:XIV  LN:784333
@SQ SN:XV   LN:1091291
@SQ SN:XVI  LN:948066
@SQ SN:Mito LN:85779

GTF:
IV  sgd gene    8683    9756    .   -   .   gene_id "YDL246C"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    8683    9756    .   -   .   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL246C_mRNA-E1";
IV  sgd CDS 8686    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL246C";
IV  sgd start_codon 9754    9756    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  8683    8685    .   -   0   gene_id "YDL246C"; transcript_id "YDL246C_mRNA"; exon_number "1"; gene_name "SOR2"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "SOR2"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    17577   18566   .   -   .   gene_id "YDL243C"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    17577   18566   .   -   .   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDL243C_mRNA-E1";
IV  sgd CDS 17580   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDL243C";
IV  sgd start_codon 18564   18566   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  17577   17579   .   -   0   gene_id "YDL243C"; transcript_id "YDL243C_mRNA"; exon_number "1"; gene_name "AAD4"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "AAD4"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd gene    1248154 1249821 .   -   .   gene_id "YDR387C"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding";
IV  sgd transcript  1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd exon    1248154 1249821 .   -   .   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; exon_id "YDR387C_mRNA-E1";
IV  sgd CDS 1248157 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding"; protein_id "YDR387C";
IV  sgd start_codon 1249819 1249821 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";
IV  sgd stop_codon  1248154 1248156 .   -   0   gene_id "YDR387C"; transcript_id "YDR387C_mRNA"; exon_number "1"; gene_name "CIN10"; gene_source "sgd"; gene_biotype "protein_coding"; transcript_name "CIN10"; transcript_source "sgd"; transcript_biotype "protein_coding";

Edit: the output of ReadsPerGene.out.tab is (thanks
swbarnes2) :
N_unmapped  73359652    73359652    73359652
N_multimapping  15789900    15789900    15789900
N_noFeature 192619  5872967 380420
N_ambiguous 648956  13479   430527

Would anyone know why most reads map to no_feature?

You are already running it with GeneCounts, which should output the same as htseq-count.  What does your ReadsPerGene.out.tab say?

Do your gtf contains exon entris (i.e. with "exon" on the 3rd column)? Could you please update your question with a few exon entries from your gtf?


text:
What do "PnBS" and "PnMS" mean in FCS flowCore parameters?
I'm using the flowCore packages to analyze cytometry data, and I am trying to understand all the parameters stored in the FCS file.
I've figured out most of the elements, but despite a thorough search I haven't been able to find documentation on some of the parameters.
For each channel, there are 9 associated parameters. Here's P1 and its associated values from one of my files :

$P1N    FSC-A
$P1R    262144
$P1B    32
$P1E    0,0
$P1V    225
$P1G    1
$P1DISPLAY  LIN
$P1BS   0
$P1MS   0

I've figured out the following from various sources (notably the FCS data standards, given here) or from the fields themselves:

$PnN Short name for parameter n.
$PnR Range for parameter number n.
$PnB Number of bits reserved for parameter number n.
$PnE Amplification type for parameter n.
$PnV Detector voltage for parameter n.
$PnG Amplifier gain used for acquisition of parameter n.
$PnDISPLAY  the transformation to apply when displaying the data (e.g. LIN, LOG)

However, these 2 don't appear in the file standard description, and flowCore does not have documentation on the parameters either.

$PnBS
$PnMS

I've seen $PnBS set at -1, 0, 48, 57, 54, 52 or 25 for various channels in my data, while $PnMS always seems to be at 0.
Does someone know what it is, or where to look for a definition?
If that helps, my FCS files were generated through BD FACSDiva Software Version 9.0.1.

As this is a package from Bioconductor I would recommend to ask on their site: support.bioconductor.org. Probably you'll get a helpful answer, but it would be nice to know where do you find these parameters

I found my answer in an old mailing list. Credit to Josef Spidlen who authored the mail.
Short answer: both parameters are proprietary keywords from the BD FACSDiva software for storing scaling values.
Long answer from the mailing list:

PnBS Indicates the automatically set Biexponential scale value for Instrument  Parameter (-1 indicates an invalid scale value). When
importing an FCS  file, if the FCS file has valid Biexponential
scales, they should be  applied to the data. If the FCS file doesn't
have valid biexponential  scales then the scales are calculated from
the entire data. If  biexponential scales exist, then if AUTOBS is
set, the PnBS values are  used and applied to the data. If AUTOBS is
not set, the PnMS values are  used and applied to the data. The
information represented by this keyword becomes part of a Tube's
attributes, and does not appear in the Tube's list of keywords. It can
be viewed in the Biexponential inspector' when the auto scaling option
is selected.
PnMS Indicates the manually set Biexponential scale value for Instrument  Parameter (-1 indicates an invalid scale value). When
importing an FCS  file, if the FCS file has valid Biexponential
scales, they should be  applied to the data. If the FCS file doesn't
have valid biexponential  scales then the scales are calculated from
the entire data. If  biexponential scales exist, then if AUTOBS is
set, the PnBS values are  used and applied to the data. If AUTOBS is
not set, the PnMS values are  used and applied to the data. The
information represented by this keyword becomes part of a Tube's
attributes, and does not appear in the Tube's list of keywords. It can
be viewed in the Biexponential inspector' when the manual scaling
option  is selected.


text:
MetaQuast for assembling samples from complex communities
I'm working with whole genome metagenomic samples from human skin, and I'm using MEGAHIT for assembly and MetaQuast for evaluation. However, MetaQuast requires a list of reference genomes for the analysis (if I understood correctly). When I perform the evaluation without the references, it uses BLASTN for aligning contigs to SILVA 16S rRNA database. Since I'm working with whole genome sequences, this is not an optimal choice. As expected, MetaQuast is identifying only two species in my contigs file, although I have manually checked through BLASTN that there are contigs that match other species' genomes. And this evaluation makes my assembly looks terrible. Since I'm working with samples from human skin, it wouldn't be viable to download a database of reference genomes and perform MetaQuast using a local alignment. So, what would be the best approach to perform the evaluation correctly in my case?
Looking forward to discussing it!

I would suggest first using a binning approach such as MetaBAT2 (or an equivalent tool, there are many) to separate your metagenomic assembly into bins roughly corresponding to genomes. I would then run CheckM to assess the bins' quality/completeness/purity.
If you still want to use quast, you could then run metaquast with a reference genome appropriate to each bin independently. I would suggest GTDB-tk to perform a classification and identify such an appropriate reference genome.

Thanks @MaximilianPress good way to take the question forward

Thank you @MaximilianPress, I received a similar answer at Reddit, so I believe this might be a better way to proceed indeed


text:
Calculate interval of flanking introns
I am looking for a easy way to calculate the intervals of flanking introns based on a transcript annotation file. I have chosen a very complicated method:
bedtools subtract -s -a genes.bed -b onlyExons.bed > introns_hg38.bed
Then I calculated the start and end positions of the transcripts from transcript.gtf annotation file and searched for adjacent intron intervals in R. Unfortunately I don't find all the flanking introns. Perhaps this is due to the 3'UTR or 5'UTR or other sequences that are not introns or exons. Any suggestions?


text:
Gene/protein expression specific to a group in omics
I am wondering what is the significance of finding a particular protein specific to a disease or control group? when we detect 1000s of proteins in a proteomics experiment, how can one be sure that the specific proteins found in say, a disease group is really absent in control group than by random chance?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.

You can never be 100% sure that a difference is real and not a statistical fluke. This is the entire concept underlying Type I errors - these are errors where the null hypothesis (usually, that there is no difference) is incorrectly rejected (a difference is called significant, when in fact none truly exists). Especially when performing many hypothesis tests, it's possible that a variable with truly identical distributions between two groups appears to have very different distributions simply due to sampling variability.
You can reduce the chance of making such an error by reducing the level of alpha used to determine significance, and by performing appropriate multiple hypothesis correction. With stringent significance thresholds and conservative correction, you can make it extremely unlikely that a "significant" difference could be due to random chance, but it's not usually possible to conduct an experiment with a Type I error rate of 0%.

But since the proteins are specific to group x , and therefore lacks any measurement in the group y (since they are not detected), Thus the scenario here is not to see if the means are different. how to run a statistical significance test against something not detected ( should it be assumed 0 mean and 0 sd for the non-detected group? and then to determine with the desired p-value threshold?)


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

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

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

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


text:
How to assess quality of WGCNA module identification in blockwisemodules()
I'm exploring WGCNA for bulkRNA sequencing analysis with human subjects.  I have healthy, myocarditis, and heart failure patients (which can be further broken into ischemic or nonischemic)
I have been spending quite a bit of time tinkering with the parameters in blockwisemodules() to identify modules/MEs, and I'm only able to assess the quality of this step with a rudimentary intuition.
For example, if I use
net = blockwiseModules(datExpr, power = 10,maxBlockSize = 20000,TOMType = "unsigned", minModuleSize = 20,reassignThreshold = 50, mergeCutHeight = 0.35,numericLabels = TRUE, pamStage = TRUE, pamRespectsDendro = TRUE)
I get what is so far the "best" module detection - many modules (although small), and about 60-70% of genes seem to fall into a module. Image below.
I find myself at a bit of a loss as to how to judge the quality of this step.  The module detection looks extremely messy compared with the tutorial data, which i wasn't surprised at really, but can anyone provide guidance here?
I can provide other images of what I get from this step that looks like much poorer detection of modules if it would help.
Last related question, I have 19 healthy, and about 41 patient samples, I wasn't sure if the best approach to module detection is to put all samples together, then look at what modules are associated with disease states which i could convert into numeric for that step, or would it be better as i have started - to detect modules in healthy, or patient groups in independent datasets?


text:
Which components must be included in the HDR plasmid?
I am trying to design a plasmid for our knock in experiment and wanted to clarify some technical issues.
I wanted to know exactly which components must be included in the plasmid.
I know that I have to include the 5' homology arm , insertion gene sequence, marker sequence and 3' homology arm.
But, I would like to know for certain the following details:

Should I take the whole sequence of the gene or do I need to remove the introns/UTR regions or any other unnecessary part of the sequence?
Which sequence should I choose for gene insertion , gene sequence or mRNA transcript sequence?
When designing plasmids, should I insert extra promoters and/or terminators with the insertion genes?

Additionally, can you recommend an online tool to design HDR Plasmid?
I have already tried using Alt R HDR design tool and Benchling but our target genome is not available there.

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

You would probably get better answers on the Biology forum. It would also help to mention whether you are using a commercial backbone for the plasmid or designing the whole construct from scratch.

Thank you @Cloudberry for your suggesion.


text:
How To display bar graph with error bars?
I am using the below code to get the bar plot. I want to get the bar plot with error bars displayed.
import seaborn as sns

def plot_coefficients(target, coeffs_df, n=10):
    if coeffs_df.shape[0]==0:
        return None # nothing to plot here!
    data = coeffs_df.copy()
    data = data.sort_values(('Coefficient','mean'), ascending=False).head(n)
    g = sns.barplot(
        data=data,
        y='Gene',
        x=('Coefficient','mean'),
        hue='Feature Type',
        hue_order=[
            h for h in [
                'Genomic Variant',
                'Gene Expression',
                'Proteomics',
                'Methylation',
            ]
        ],
        dodge=False
    )
    g.set(
        xlim=(0,1),
        xlabel='Importance [0,1]',
        ylabel='',
        title=f"Top {target} predictor features"
    )
    return g

Coefficient here is a score generated from predictive modelling that ranges from 0 to 1. 2. data here is entire data and I created a mean column. Ideally I wanted to calculate it in seaborn plot rather than creating a mean column separately before.

I presume you import seaborn as sns, right?

@KamilSJaron Yes

I did put forward a comment, it was removed I think the idea is its a full post or zero. The problem proving a negative is tricky.
I don't think its doable within the exact system specified, i.e. a Seaborn bar chart. Again I could be wrong.
If the OP went to a Seaborn dot line formation - that would work sns.barplot is changed to sns.pointplot but then all the nice colouring disappears. Standard deviation  is one error range if all the values are present: but I think the OP will need to request it (.sd() ?? can't remember).
In matplotlib there is an error bar method where the coder specifies the error bars and that is plt.errorbar. An example from the official docs is here. Thus call a bar chart and then add the error bars onto it.
Summary The reasons I don't think its doable in Seaborn:

I'm pretty certain Seaborn has no independent errorbar method unlike matplotlib
Seaborn can do error bars given the full data set, but those are created automatically for a given method. Maybe there is a "sns.barplot-errorbar" method but I don't know it.

The Seaborn approach does carry the advantage of not having to work out the error bar independently.


text:
snpEff annotation and variant in a gene
I'm using snpEff to annotate genetic variants. However, the annotation often returns more than one annotation (see below for an example):
NZ_JAKVDJ010000001.1    119361  .   T   C   23.1481 .   DP=22;VDB=0.124119;SGB=-0.690438;RPB=0.369955;MQB=0.0631559;MQSB=0.000399184;BQB=1;MQ0F=0;AC=1;AN=1;DP4=0,5,12,5;MQ=5;ANN=C|synonymous_variant|LOW|oprB|gene-MMZ02_RS00620|transcript|TRANSCRIPT_gene-MMZ02_RS00620|protein_coding|1/1|c.444A>G|p.Leu148Leu|444/1359|444/1359|148/452||WARNING_TRANSCRIPT_NO_START_CODON,C|upstream_gene_variant|MODIFIER|MMZ02_RS00610|gene-MMZ02_RS00610|transcript|TRANSCRIPT_gene-MMZ02_RS00610|protein_coding||c.-3612A>G|||||3612|,C|upstream_gene_variant|MODIFIER|MMZ02_RS00615|gene-MMZ02_RS00615|transcript|TRANSCRIPT_gene-MMZ02_RS00615|protein_coding||c.-1069A>G|||||1069|,C|downstream_gene_variant|MODIFIER|MMZ02_RS00625|gene-MMZ02_RS00625|transcript|TRANSCRIPT_gene-MMZ02_RS00625|protein_coding||c.*785A>G|||||785|,C|downstream_gene_variant|MODIFIER|MMZ02_RS00630|gene-MMZ02_RS00630|transcript|TRANSCRIPT_gene-MMZ02_RS00630|protein_coding||c.*1063A>G|||||1063|,C|downstream_gene_variant|MODIFIER|MMZ02_RS00635|gene-MMZ02_RS00635|transcript|TRANSCRIPT_gene-MMZ02_RS00635|protein_coding||c.*2025A>G|||||2025|,C|downstream_gene_variant|MODIFIER|MMZ02_RS00640|gene-MMZ02_RS00640|transcript|TRANSCRIPT_gene-MMZ02_RS00640|protein_coding||c.*2876A>G|||||2876|,C|downstream_gene_variant|MODIFIER|MMZ02_RS00645|gene-MMZ02_RS00645|transcript|TRANSCRIPT_gene-MMZ02_RS00645|protein_coding||c.*3698A>G|||||3698|,C|intragenic_variant|MODIFIER|gene-MMZ02_RS01495|null|gene_variant|null|||n.119361T>C||||||    GT:PL   1:50,0

However, now I wonder if the first annotation always refers to the observed variant (whose position is listed at the start of the line with the annotations an other quality parameters. Thus, if I got the annotation from the example above, can I always trust that the variant I'm looking at occurred in the gene oprB? Or will snpEff sometimes also give annotations where only the 2nd, 3rd, etc. annotation is directly related to the varying DNA position? If so, how can I know with certainty in which gene one variant appears?
Thanks a lot for your help in advance!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.


text:
Bacterial genomes and snpeff warnings 'WARNING_TRANSCRIPT_NO_START_CODON'
I annotated a bunch of genomes with snpeff. And although the original reads' qualities and the matches with the reference genomes were quite high, I see the above mentioned warning in almost every variant in my vcf file. Thus, I wonder how this could happen. Since I'm using bacterial genomes with a manually built database (from these reference genomes) and there are quite some genes (especially in gene clusters) that have the alternative start codon GTG, I suspect this causing the warning. Any idea how to test if this is really the reason for the warning?
Thanks a lot in advance!
Edit:
As requested, here is the html snpEff report.
And I've also uploaded an example vcf file after annotation here

Bacteria have a different translation table from eukaryote nuclear DNA; it may be that SnpEff is expecting human genes in the genome. Have you followed the build steps outlined in the documentation?
https://pcingola.github.io/SnpEff/se_build_db/

If your genome, or a chromosome, uses non-standard codon tables you must update snpEff.config accordingly. A typical case is when you use mitochondrial DNA. Then you specify that chromosome 'MT' uses codon. Invertebrate_Mitochondrial codon table. Another common case is when you are adding a bacterial genome, then you specify that the codon table is Bacterial_and_Plant_Plastid.

Thanks for your suggestion. You're of course right, I should have done this before. And now the WARNING_TRANSCRIPT_NO_START_CODONs disappeared, but other warnings appeared ... (WARNING_MULTIPLE_STOP_CODONS & WARNING_TRANSCRIPT_INCOMPLETE).

According to your report, you are using human HG18 as the reference genome version.  Therefore, the report will be wrong along with the wrong codons warnings.  It is unfruitful to try with a human database with bacteria.

Yes, thanks for your comment - already changed it yesterday. It was indeed the problem.


text:
Handling of sex chromosomes in male when computing CNV copy ratio
Reading from CNVkit documentation, Mathworks documentation, Bionano documentation, for CNV calling, it seems that the common practice for array is to compute the linear copy ratio by normalizing against the expected depth for the reference ploidy, i.e.,:
$$CR = \frac{D_{OBS}}{P_{REF}}$$
where $D_{OBS}$ is the observed depth, and $P_{REF}$ is expected depth for the reference ploidy for a specific segment (i.e., expected depth for 2 copies for females, and males outside sex chromosomes or in PAR regions).
Is that correct? Are there alternative approaches used in array copy number calling?


text:
Count reads at specific gene features
I have BAM files from an RNA-Seq experiment and for all genes (or a subset) I want to get the number of reads in regions around the TSS (e.g. 2kbp) and the TES (e.g. 2 kbp) and calculate the ratio between them. I was thinking of first creating these regions for each gene and store them in a BED file and then use bedtools intersect with the BAM files to get the reads in these regions. My questions is how would I make such a BED file that contains only the specified regions downstream of the TSS and specified regions upstream of the TSS?

First of all, you need gene coordinates. You can retrieve them in GTF format from GENCODE or other sources.
Next, extract TSS and TES +/- 2kb windows from each gene, taking in mind that genes can be in the forward or reverse strand. Something like this:
awk -F"\t|\"" -vOFS="\t" '$3=="gene" {
  TSS=($7=="+" ? $4-2000"\t"$4+2000 : $5-2000"\t"$5+2000)
  TES=($7=="+" ? $5-2000"\t"$5+2000 : $4-2000"\t"$4+2000)
  print $1,TSS,$10"_TSS_2kb",0,$7
  print $1,TES,$10"_TES_2kb",0,$7
}' annotation.gtf > tss_tes.bed

Here I used GNU Awk from Bash, but you can write your own solution using Python, R or whatever you like.
Finally, count reads falling in these intervals with your tool of choice (you mentioned bedtools, which is fine).


text:
Find a genomic coordinates for a protein aminoacid position
Is there a function that can map genomic position (hg19) back from a protein position? I can have name of a particular transcript, and exon number.
For example, I have KRAS gene for which I would like to know the exact genomic position of G12, and I expect to have an output saying:
chr12, from 25398283 to 25398285.
transvar could be an excellent solution if the current built was working (broken links for fasta annotations). The web search is working but not allowing files heavier than 5kB which is my case.

This is considerably more complicated since there can be more than one genomic position referring to G12 in the protein in cases where multiple transcripts use different exons and yet >1 of them result in a protein with a G at position 12.

@lizaveta could you please review the answer? In my opinion this does answer your question (just my opinion), if it does please kindly consider "accepting" the answer.

You can do this manually by looking up the protein variant on https://varsome.com. For your example, that would mean searching for KRAS G12V (change V for anything else, the idea is to give an HGVS protein notation for a variant, so any single letter amino acid code would work). That would take you to this page https://varsome.com/variant/hg19/KRAS%20G12V?annotation-mode=germline, which includes:

Of course, I imagine you want to automate this. VarSome can help there too, through our API, but that is a paid service. If you have the option of paying for this, the VarSome API can be used to (among other things) look up arbitrary variants and translate between different notations, so you could easily adapt this to get the information you need.
Disclosure: I work for the company behind VarSome.


text:
Parsing text within a taxonomy.csv file
I have the following file in tab separated values (TSV).

ID
Taxonomy

76980f6d906c6baaef3f96559fc7ba1b
d__Eukaryota; p__Cercozoa; c__Vampyrellidae; o__Vampyrellidae; f__Vampyrellidae; g__uncultured

4f808d3adec894395304661391037369
d__Eukaryota; p__Ochrophyta; c__Xanthophyceae; o__uncultured; f__uncultured; g__uncultured

777584465097568f5b30bc9e2b6584c1
d__Eukaryota; p__Labyrinthulomycetes; c__uncultured; o__uncultured; f__uncultured; g__uncultured

5a99dafabe78847d8f5fed0d37f8672a
d__Eukaryota; p__Labyrinthulomycetes; c__Labyrinthulomycetes; o__Labyrinthulomycetes; f__Thraustochytriaceae; g__uncultured

I would like to modify the "g__uncultured" with the last known rank possibly using either shell or R scripts.
The new table would then read:

ID
Taxonomy

76980f6d906c6baaef3f96559fc7ba1b
d__Eukaryota; p__Cercozoa; c__Vampyrellidae; o__Vampyrellidae; f__Vampyrellidae; g__f__Vampyrellidae

4f808d3adec894395304661391037369
d__Eukaryota; p__Ochrophyta; c__Xanthophyceae; o__uncultured; f__uncultured; g__c__Xanthophyceae

777584465097568f5b30bc9e2b6584c1
d__Eukaryota; p__Labyrinthulomycetes; c__uncultured; o__uncultured; f__uncultured; g__p__Labyrinthulomycetes

5a99dafabe78847d8f5fed0d37f8672a
d__Eukaryota; p__Labyrinthulomycetes; c__Labyrinthulomycetes; o__Labyrinthulomycetes; f__Thraustochytriaceae; g__f__Thraustochytriaceae

It would help if you could provide your attempts to solve this, it would be easier to help you fix the issue than propose an answer.

The pandas solution is here:
#!/usr/bin/env python3
import pandas as pd
from pathlib import Path

path = '/Whatever'
infile = 'myfile.csv'
outfile = 'outfile.csv'
pd.set_option('display.max_colwidth', None)
df = pd.read_csv(Path(path,infile), sep='\t')
taxaChange = {'76980f6d906c6baaef3f96559fc7ba1b':'g__f__Vampyrellidae', '4f808d3adec894395304661391037369':'g__c__Xanthophyceae'}
df2 = df.set_index('ID')
df3 = df2
for x,y in taxaChange.items():
    df3.loc[x, 'Taxonomy'] = df2.loc[x, 'Taxonomy'].replace('g__uncultured', y)
print (df3)

df.to_csv(Path(path,outfile), sep='\t')

Output
ID                                       Taxonomy
76980f6d906c6baaef3f96559fc7ba1b         d__Eukaryota; p__Cercozoa; c__Vampyrellidae; o__Vampyrellidae; f__Vampyrellidae; g__f__Vampyrellidae
4f808d3adec894395304661391037369         d__Eukaryota; p__Ochrophyta; c__Xanthophyceae; o__uncultured; f__uncultured; g__c__Xanthophyceae
777584465097568f5b30bc9e2b6584c1         d__Eukaryota; p__Labyrinthulomycetes; c__uncultured; o__uncultured; f__uncultured; g__uncultured
5a99dafabe78847d8f5fed0d37f8672a         d__Eukaryota; p__Labyrinthulomycetes; c__Labyrinthulomycetes; o__Labyrinthulomycetes; f__Thraustochytriaceae; g__uncultured

Both g__f__Vampyrellidae and g__c__Xanthophyceae are there, i.e. the code was tested and worked fine.
The point for the pandas approach is you've a natural table format so it makes sense to take it straight to a dataframe (table -> dataframe). regex-ing using \t and \n is of course the alternative and you need to formulate the table as a hash (Perl) or dictionary (Python).
Notes The df and df2 can reduced a single dataframe by updating the df (inplace was omitted and will work for an index), rather than creating a new df. However, the df3 is needed (need to copy the dataframe) because 'inplace' didn't work. inplace is a command to force the dataframe to update. The set_option is just so I can print the display without truncation.
Finally
You substitute the infile and outfile to correct path and name values, save it as a .py file and just python3 myscript.py.

Thanks for your reply. However, my tsv file was only an example. I've a lot of g__uncultured that should be replaced. Hence, the taxaChange should be something more general, able to parse different lines.

Hi @Marco. This does answer the specific question in my view and the code was tested and worked fine. It is 70% towards a generalised solution. However, automating the code to process lots of stuff is a separate question as general principle and particularly because its not clear how you want to do that (there are no indicators except for the above comment). I suspect you've got something is specific in mind. You can use this code as the basis for that question if you wish to pursue a pandas solution. pandas is "trendy", but there are lots of ways to do this.


text:
RNA strings and amino acid string
Which of the following RNA strings could translate into the amino acid string PRTEIN? (Select all that apply.)
I am not sure if the following mass of amino acid is relevant to the question.

CCAAGUACAGAGAUUAAC

CCAAGAACAGAUAUCAAU

CCACGUACUGAAAUUAAC

CCGAGGACCGAAAUCAAC

I don't quite get the question but here's the answer and the code
#!/usr/bin/env python3
from Bio.Seq import Seq

for dna in ['CCAAGUACAGAGAUUAAC','CCAAGAACAGAUAUCAAU', 'CCACGUACUGAAAUUAAC','CCGAGGACCGAAAUCAAC']:
    print(Seq(dna).translate())

PSTEIN PRTDIN PRTEIN PRTEIN

may I know where is the API code documentation for .translate() function call ?

Its Biopython part of the Bio package. Its not an API its a Python package

https://biopython.org/docs/1.75/api/Bio.Seq.html#Bio.Seq.Seq.translate

what is the opposite of .translate() function calls ?  I mean let's say I am given an amino acid string CYCLIC, how do I obtain all the possible combination of DNA strings before the translation process ?

Thats a separate question because it will involve degeneracy.

ok, let me ask in a separate question thread

https://bioinformatics.stackexchange.com/questions/19498/reverse-translation-from-amino-acid-string-to-dna-strings

Minor note: A sequence containing U usually means it is an RNA sequence, while T would hint at it being a DNA sequence. In eukaryotes, the DNA sequence is often not directly translated into the corresponding RNA sequence, because of effects like alternative splicing.

Thanks @RobertHass. Could be a virus?

Here's how you can do the translation yourself, using the genetic code in form of a dictionary from Biopython:
import Bio

# Load the genetic code (a mapping from RNA triplet to amino acid)
d = Bio.Data.CodonTable.standard_rna_table.forward_table
print('Genetic code:', d)
print()

# Use it to translate an RNA sequence into an amino acid sequence
s = 'CCAAGUACAGAGAUUAAC'
for triplet in (s[i:i+3] for i in range(0, len(s), 3)):
    print(triplet, '->', d[triplet])

Printed result:
Genetic code: {'UUU': 'F', 'UUC': 'F', 'UUA': 'L', 'UUG': 'L', 'UCU': 'S', 'UCC': 'S', 'UCA': 'S', 'UCG': 'S', 'UAU': 'Y', 'UAC': 'Y', 'UGU': 'C', 'UGC': 'C', 'UGG': 'W', 'CUU': 'L', 'CUC': 'L', 'CUA': 'L', 'CUG': 'L', 'CCU': 'P', 'CCC': 'P', 'CCA': 'P', 'CCG': 'P', 'CAU': 'H', 'CAC': 'H', 'CAA': 'Q', 'CAG': 'Q', 'CGU': 'R', 'CGC': 'R', 'CGA': 'R', 'CGG': 'R', 'AUU': 'I', 'AUC': 'I', 'AUA': 'I', 'AUG': 'M', 'ACU': 'T', 'ACC': 'T', 'ACA': 'T', 'ACG': 'T', 'AAU': 'N', 'AAC': 'N', 'AAA': 'K', 'AAG': 'K', 'AGU': 'S', 'AGC': 'S', 'AGA': 'R', 'AGG': 'R', 'GUU': 'V', 'GUC': 'V', 'GUA': 'V', 'GUG': 'V', 'GCU': 'A', 'GCC': 'A', 'GCA': 'A', 'GCG': 'A', 'GAU': 'D', 'GAC': 'D', 'GAA': 'E', 'GAG': 'E', 'GGU': 'G', 'GGC': 'G', 'GGA': 'G', 'GGG': 'G'}

CCA -> P
AGU -> S
ACA -> T
GAG -> E
AUU -> I
AAC -> N

Note that the same RNA sequence can also be translated into another sequence of amino acids by shifting the reading frame 1 or 2 positions. For example, the RNA sequence 'CCAAGUACAGAGAUUAAC' can be split into different triplets:

'CCA', 'AGU', 'ACA', ...
'CAA', 'GUA', 'CAG', ...
'AAG', 'UAC', 'AGA', ...

In biological systems, an mRNA would have a start codon that marks the beginning of an open reading frame, which is often the triplet AUG.


text:
reverse translation from amino acid string to DNA strings
what is the opposite of .translate() function calls ?
I mean let's say I am given an amino acid string CYCLIC, how do I obtain all the possible combinations of DNA strings before the translation process ?
Note: I was told that degeneracy has an important role in this particular aspect, but I am not entirely sure about the degeneracy calculation in the following screenshot highlighted in yellow color.

The degenerate code is as follows:
table = {
        'ATA':'I', 'ATC':'I', 'ATT':'I', 'ATG':'M',
        'ACA':'T', 'ACC':'T', 'ACG':'T', 'ACT':'T',
        'AAC':'N', 'AAT':'N', 'AAA':'K', 'AAG':'K',
        'AGC':'S', 'AGT':'S', 'AGA':'R', 'AGG':'R',
        'CTA':'L', 'CTC':'L', 'CTG':'L', 'CTT':'L',
        'CCA':'P', 'CCC':'P', 'CCG':'P', 'CCT':'P',
        'CAC':'H', 'CAT':'H', 'CAA':'Q', 'CAG':'Q',
        'CGA':'R', 'CGC':'R', 'CGG':'R', 'CGT':'R',
        'GTA':'V', 'GTC':'V', 'GTG':'V', 'GTT':'V',
        'GCA':'A', 'GCC':'A', 'GCG':'A', 'GCT':'A',
        'GAC':'D', 'GAT':'D', 'GAA':'E', 'GAG':'E',
        'GGA':'G', 'GGC':'G', 'GGG':'G', 'GGT':'G',
        'TCA':'S', 'TCC':'S', 'TCG':'S', 'TCT':'S',
        'TTC':'F', 'TTT':'F', 'TTA':'L', 'TTG':'L',
        'TAC':'Y', 'TAT':'Y', 'TAA':'_', 'TAG':'_',
        'TGC':'C', 'TGT':'C', 'TGA':'_', 'TGG':'W',
    }

There isn't a function that will do the reverse translation, so it user code and requires the degenerate amino acid table. Degeneracy is important into primer design agreed. I'm not around for several days.

May I know where exactly did you obtain the degenerate code table ?

There everywhere, I got it from Geek for Geeks

Given the table above the solution for CYCLE is:
from itertools import product

prot = 'CYCLE'
degenDict = dict()
for v,k in table.items():
    if k in degenDict:
        degenDict[k].append(v)
    else:
        degenDict[k]=[]
        degenDict[k].append(v)
nucs = [degenDict[resid] for resid in prot]
for degenNuc in product(*nucs):
    print(''.join(degenNuc))

TGCTACTGCCTAGAA
TGCTACTGCCTAGAG
TGCTACTGCCTCGAA
TGCTACTGCCTCGAG
TGCTACTGCCTGGAA
TGCTACTGCCTGGAG
TGCTACTGCCTTGAA
TGCTACTGCCTTGAG
TGCTACTGCTTAGAA
TGCTACTGCTTAGAG
TGCTACTGCTTGGAA
TGCTACTGCTTGGAG
TGCTACTGTCTAGAA
TGCTACTGTCTAGAG
TGCTACTGTCTCGAA
TGCTACTGTCTCGAG
TGCTACTGTCTGGAA
TGCTACTGTCTGGAG
TGCTACTGTCTTGAA
TGCTACTGTCTTGAG
TGCTACTGTTTAGAA
TGCTACTGTTTAGAG
TGCTACTGTTTGGAA
TGCTACTGTTTGGAG
TGCTATTGCCTAGAA
TGCTATTGCCTAGAG
TGCTATTGCCTCGAA
TGCTATTGCCTCGAG
TGCTATTGCCTGGAA
TGCTATTGCCTGGAG
TGCTATTGCCTTGAA
TGCTATTGCCTTGAG
TGCTATTGCTTAGAA
TGCTATTGCTTAGAG
TGCTATTGCTTGGAA
TGCTATTGCTTGGAG
TGCTATTGTCTAGAA
TGCTATTGTCTAGAG
TGCTATTGTCTCGAA
TGCTATTGTCTCGAG
TGCTATTGTCTGGAA
TGCTATTGTCTGGAG
TGCTATTGTCTTGAA
TGCTATTGTCTTGAG
TGCTATTGTTTAGAA
TGCTATTGTTTAGAG
TGCTATTGTTTGGAA
TGCTATTGTTTGGAG
TGTTACTGCCTAGAA
TGTTACTGCCTAGAG
TGTTACTGCCTCGAA
TGTTACTGCCTCGAG
TGTTACTGCCTGGAA
TGTTACTGCCTGGAG
TGTTACTGCCTTGAA
TGTTACTGCCTTGAG
TGTTACTGCTTAGAA
TGTTACTGCTTAGAG
TGTTACTGCTTGGAA
TGTTACTGCTTGGAG
TGTTACTGTCTAGAA
TGTTACTGTCTAGAG
TGTTACTGTCTCGAA
TGTTACTGTCTCGAG
TGTTACTGTCTGGAA
TGTTACTGTCTGGAG
TGTTACTGTCTTGAA
TGTTACTGTCTTGAG
TGTTACTGTTTAGAA
TGTTACTGTTTAGAG
TGTTACTGTTTGGAA
TGTTACTGTTTGGAG
TGTTATTGCCTAGAA
TGTTATTGCCTAGAG
TGTTATTGCCTCGAA
TGTTATTGCCTCGAG
TGTTATTGCCTGGAA
TGTTATTGCCTGGAG
TGTTATTGCCTTGAA
TGTTATTGCCTTGAG
TGTTATTGCTTAGAA
TGTTATTGCTTAGAG
TGTTATTGCTTGGAA
TGTTATTGCTTGGAG
TGTTATTGTCTAGAA
TGTTATTGTCTAGAG
TGTTATTGTCTCGAA
TGTTATTGTCTCGAG
TGTTATTGTCTGGAA
TGTTATTGTCTGGAG
TGTTATTGTCTTGAA
TGTTATTGTCTTGAG
TGTTATTGTTTAGAA
TGTTATTGTTTAGAG
TGTTATTGTTTGGAA
TGTTATTGTTTGGAG


text:
mass spectrum matching with cyclic peptide and linear peptide
I was doing a simple quiz which I encountered the following two questions related to mass spectrum matching with cyclic peptide and linear peptide.
I have done some preliminary google search, but I am not sure how to use either  https://codereview.stackexchange.com/questions/193859/find-a-cyclic-peptide-with-a-given-non-ideal-spectrum or https://github.com/rochoa85/pepfun in this case.

Which of the following cyclic peptides could have generated the
theoretical spectrum 0 71 101 113 131 184 202 214 232 285 303 315 345
416? (Select all that apply.)

TAIM

TMIA

MAIT

MTAI

TMLA

TLAM

Which of the following linear peptides is consistent with Spectrum =
{0 71 99 101 103 128 129 199 200 204 227 230 231 298 303 328 330 332
333}? (Select all that apply.)

AVQ

TCQ

TVQ

QCV

TCE

ETC

You are way overthinking this.
As this is homework I won't give any code, but I will throw back some questions:

What is the definition of m/z?
Were the only tool available Wikipedia, could you simply read the infobox of each residue?
If so, what values would you need and what elements do you need to subtract?
Does residue order matter for linear and backbone cyclised peptides? (Let's ignore isoaspartate etc.)
Have you been introduced to RDKit in your course? E.g. what would this give you?

from rdkit import Chem
from rdkit.Chem import Descriptors
mol: Chem.Mol = Chem.MolFromFASTA('ELVIS')
print( Descriptors.ExactMolWt(mol) )

Sorry for late reply.  I am stucked with many other things last week.  I am new to rdkit, but  print( Descriptors.ExactMolWt(mol) ) would give 559.32172802 , and what can I learn from this number ?

and what do you exactly mean by infobox of each residue ?

https://en.wikipedia.org/wiki/Mass-to-charge_ratio and https://en.wikipedia.org/wiki/Help:Infobox

May I ask why you specifically chose ELVIS ?  or was it just a demo example ?

and I am not sure about Does residue order matter for linear and backbone cyclised peptides? (Let's ignore isoaspartate etc.)

Demo. The amino acid sequence ELVIS is an ancient joke as it is very common, but ELVISISALIVE was not present in the draft human genome.

but there is difference between linear peptide and cyclic peptide (which is more difficult to solve) ?

Wait, may I ask why would I need https://en.wikipedia.org/wiki/Mass-to-charge_ratio for this question ?


text:
Where does the choice of k = 51 for de Bruijn graphs of large genomes come from?
This question has also been asked on Biostars
From hearsay I know that de Bruijn graphs of large genomes (e.g. human) are usually constructed with k = 51, or that k = 51 is at least a good initial choice.
I however am unable to find any source for this, does anyone know where it is coming from?

Not an answer. SGA is the first assembler I know that has a 51-mer default (more exactly ≥51bp overlap) but I am not sure if k=51 comes from SGA. Note that ALLPATHS-LG/DISCOVAR, arguably the best short-read dBG assembler for large genomes, uses k>100 for end-merged reads, so I don't think k=51 is necessarily an optimal choice.

Note that k should (must?) be odd to avoid kmers that are reverse-complement palindromes.