text:
Are variant calling files personally identifiable information?
This question is to be read in the context of data protection. It seems to be common understanding that the whole genome sequence of an individual is personally identifiable, non-anonymizable information: if you possess the whole genome sequence, you can without any doubt identify an individual. You also cannot anonymize it, because you would need to obstruct the genome sequence in some way, making it useless for any diagnostics or research.
But how about Variant Calling files (specifically referring to the VCF file format)? You reference against a known genome, but you neglect structural changes such as translocations, duplications and deletions. Can you still identify a person based on a Variant Calling file?

What makes you think a VCF cannot contain structural variants? Of course it can. This doesn't really affect the answer to your question, a VCF of small vairants from even a whole exome analysis is more than enough to identify an individual under certain circumstances, but there is no reason why the VCF cannot also contain translocations and indels.

@terdon My superficial research. If a VCF contains also structural variants, the answer to my question is even clearer.

TLDR: yes! be careful with someone's genomic data!
There are two aspects to this question: can I find a match for a random VCF in a database of genomes (YES) or can I identify a subject who is not in a (public) database.
But even when it may not possible to identify the subject itself, using DTC genetic testing sites may enable you to identify (close or distant) relatives (e.g. read about the Golden State Killer). Identifying relatives brings you close to identifying the individual, and also those relatives may not be pleased with learning that the variants of someone with who they share (a small or larger part) of their genome are accessible by anyone.
We are further ahead for some traits than others, but for more and more human characteristics we will be able to predict e.g. eye color, height based on genome data, or at least estimate to some degree. Genetic ancestry and biological sex are trivial. Often genomic data is from a patient with a (rare) disorder, further narrowing down the options.
In addition, there is a difference between "can someone be identified today" (which is definitely possible but maybe tough) and "given that the variants are out there (online) forever can someone be identified within 10 years" and the answer to the latter part is a loud yes.

Yes, definitely identifiable. The combination of ~80 unlinked common SNPs can be fairly unique in the entire human population, let alone the whole VCF file.
EDIT: 30 in the original answer is an underestimate. We need ~80 unlinked common SNPs to uniquely identify every individual to a low false positive rate. Here is the derivation based on the approximate solution to the birthday paradox.
For simplicity, we assume each person has a haploid genome. Then the genotype at each SNP can only be 0 or 1. Suppose there are $n$ people in the world and we have $m$ sites. On the condition that $n^2\ll 2^{m+1}$, the probability of seeing any two people having the same binary genotypes across all $m$ sites is
$$p(n,m)\approx 1-e^{-n^2/2^{m+1}}\approx n^2/2^{m+1}$$
There are $n\approx2^{33}$ people. Then $p(2^{33},72)<1\%$.
In reality, we are diploid, so each genotype has three states instead of two. This would lead to a smaller $m$ at the same false positive rate. On the other hand, selected SNPs are not totally random. This would lead to a larger $m$. It is safer to choose a larger $m$.

Not only are they theoretically identifiable based on the combination of SNPs as others have said (you are more or less identifiable from ~20 well chosen microsats) - but it's possible to be actually identified from your data using a long-range familial search, based on genetics and a bit of sleuthing.
https://science.sciencemag.org/content/362/6415/690.abstract

There should be someone at your workplace you can ask about this, but I would say, yes, since there can be SNPs in there that are private, or even a combination of public SNPs might be identifying, it should be treated as potentially identifying.


text:
Impact of merging ChIP-seq runs of the same sample on PCR duplicates identification?
I'd like to a follow up question to this question related to merging fastq files for ChIP-seq.
Let's assume we have one sample library that is re-sequenced two or three times in order to achieve a desired read coverage. The resulting output of the sequencing is 2 or 3 fastq files for one individual sample.
If one has to mark duplicates (for example using Picard's MarkDuplicates) should the sub-samples be merged at the fastq level or at the bam file level (post alignment) after flagging duplicates before the merge?
If one takes the 2-3 sub-fastq files align them and then mark the duplicates in the bam files if the sequencing is shallow one could end up with not calling reads that are actually PCR duplicates of the library preparation. While in the other scenario, a fastq merge might result in reads to be marked as "duplicates" just because they appear 2-3 times, but maybe that read is a single reads that I sequenced 2-3 times, individually in each run.
My point is that I would expect this 2 different approached to lead to different outcomes. Probably this is a negligible effect in a overall analysis, but maybe it might actually have an impact in certain biological conditions (I'm thinking of heterochromatin silencing of repeat elements in the genome). I'm going to try it out as soon as I have a change, but I was wondering if this was discussed already somewhere else (I couldn't find much relevant) or if someone has some ideas on how to tackle this.
Let's assume that the reads are 50bp single end.

If you are sequencing same library multiple times (as you state in the question), any read that appears in  multiple sequencing runs is likely to have been generated via PCR duplication from the same original biological molecule, and therefore needs to be removed. Indeed, the chance of the read being a duplicate is the same if it was sequenced twice in the same lane, or it was sequenced in two different lanes, as long as it comes from the same library.
Thus, in this situation, you should either:

Merge the fastqs, map and MarkDuplicates or
Map the separate fastqs, merge the BAM files, and then MarkDuplicates.

You should definitely merge before doing either MarkDuplicates or peak calling.

Thanks for your reply, it makes sense. However do you think is possible to quantify this "likely to have been generated via PCR duplication". Like how much more likely?

It isn't possible, not without specially designed experiments. But the important point is that irrespective of the likelihood that a duplicate is a "biological" duplicate or a "PCR" duplicate, that probability is the same between sequencing runs as it is within a sequencing run. If it doesn't make sense to deduplicate the merged data, it wouldn't make sense to deduplicate the individual files either.

source: Stack Exchange

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

sort: votes


text:
Removing cells zero for a gene from a scRNA-seq data
I have a big single-cell RNA seq data
> dput(head(new.dat[,1:10]))
structure(list(cell1 = c(0.793763840992639, 0, 1.96843530982957,
0.461736429639991, 0.717968540649498, 0), cell2 = c(3.61741696702738,
0.231662370550224, 0, 0, 0, 0), cell3 = c(4.14348883366621, 0.118161316317251,
0.08074552209482, 2.27968429766934, 0.0470313356296409, 0), cell4 = c(1.34783143327084,
0.0094666040612932, 1.14392942941128, 0.652535826921119, 0.357542816432864,
0.149587369334621), cell5 = c(1.27104023273899, 1.55185229643731,
0, 0, 0, 0.0117525723115277), cell6 = c(1.92307653575663, 0,
0, 0.319156642478379, 0, 0), cell7 = c(3.9343015424917, 0.132824589520901,
0.119679885703561, 0.772516422897241, 0.0236884909844904, 0),
    cell8 = c(3.74969491678643, 0.103404975609384, 0.0354753982873036,
    0, 0, 0), cell9 = c(1.19084857532713, 3.9213265721495, 0,
    0.0341973245272891, 0.0419122921627454, 0), cell10 = c(4.1224255501566,
    0.301871669274068, 0.0633536200981225, 0.389959552469879,
    0, 0.0405296102106492)), row.names = c("PTPRC", "MHC-II",
"ITGAM", "Ly6C", "Ly6G", "EMR1"), class = "data.frame")
>

> dim(new.dat)
[1]     33 263086
>

How I remove every columns which are zero for one gene, let's say PTPRC?

WhichCells(seurat_object made by this matrix, slot = 'counts', expression = PTPRC > 0 )

By Seurat R package this is possible

Please accept one answer if either of them solved your problem. Thank you.

Hi there, it would be nice to write more details and an actual working code.

This is basic R subsetting operation. Most subset operations give you a subset of rows based on a column, while you're doing the opposite here - but the logic is the same.
This gives you a subset of rows when the expression is based on a column
df[logical_expression_that_subsets_rows,]

Why not try this - it should give you a subset of columns based on a row value. (Untested, but you should be able to build on this):
df[, logical_expression_that_subsets_columns]

The logical expression for the first could be df$cell1 < 0.5 (or the equivalent df[,"cell1"] < 0.5. You should be able to get a logical expression for the second based on this.

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


text:
How do I access sequence and MAPQ from bam using HTSlib in C++?
I am using the "htslib/sam.h" header in C++. I need to access the SEQ field of bam reads and store each sequence in a vector of strings. I also need to access the MAPQ field of each read to check if the read passes my threshold quality of 10 (which I have stored as an integer). Does anyone know which htslib/sam.h functions can accomplish these tasks?

The MAPQ is stored in b->core.qual, given a bam1_t *b. You can access the sequence with bam_seqi(bam_get_seq(b), some_offset)(see https://github.com/samtools/htslib/blob/develop/htslib/sam.h#L317). This will get you a base at a time, so you'd need to iterate over it (the length of the sequence being in b->core.l_qseq).

Thank you! This worked


text:
Seurat clustering Methods-resolution parameter explanation
I am learning the Seurat algorithms to cluster the scRNA-seq datasets. I found this explanation, but am confused. Can someone explain it to me, "The FindClusters function implements the procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.6-1.2 typically returns good results for single cell datasets of around 3K cells". My questions are:

What is "granularity" of the downstream analysis?
How does the parameter "resolution" work?
How to know the optimum value for parameter "resolution"?

I found this paper about the resolution that might help, but not clear to me yet. https://iopscience.iop.org/article/10.1088/1742-5468/2008/10/P10008/pdf.
I think the higher the resolution, the more clusters produced meaning the smaller clusters obtained? The smaller the clusters the more specific cell types we get, is that correct?

What is "granularity" of the downstream analysis?

Low granularity implies basic cell types (e.g. T cells, B cells), implies few larger clusters.
High granularity implies more refined cell types (e.g. CD4+ Central Memory T, Treg, Memory B cell etc.), implies more smaller clusters

How does the parameter "resolution" work?

A high resolution parameter value produces more smaller clusters, i.e. higher granularity. If you want to know how modularity clustering works, I recommend the Leiden algorithm paper https://www.nature.com/articles/s41598-019-41695-z

How to know the optimum value for parameter "resolution"?

Tricky, and it depends on what questions you want to answer. Normally it's a subjective iterative process of trying a clustering resolution, checking if the clusters make sense, re-trying etc. but I also quite like the clustree approach https://cran.r-project.org/web/packages/clustree/index.html

The smaller the clusters the more specific cell types we get, is that correct?

Yes correct


text:
Can HTSlib extract bam reads occurring in a specific region without iterating through the whole file?
I am using htslib/sam.h to write a C++ program. As part of this program, I must extract reads occurring on specific scaffolds in specific regions from a bam file. Essentially, I want to perform the following samtools command in C++:
samtools view name.bam -o name.sam scaffold:pos-pos

Since I must extract reads from thousands of regions, I do not want to iterate through the whole bam file each time I extract reads from a region. I am wondering if HTSlib can extract reads occurring in a specified scaffold and range without iterating through the whole bam. I found the function bam_itr_querys() in the htslib/sam.h documentation on GitHub, but don't know how to format the inputs for this function, or what the "idx" input is. I also saw somewhere that I would need to pair this function with bam_itr_next(), but don't know the format to do so.
EDIT: My bam files are sorted and indexed, if that matters.

Have you looked at how samtools view works with a bed file of coordinates as input?

@swbarnes2 I have not, but my issue is that I cannot use 'samtools view'. I should provide some context (I will edit the original question to include this info): in my first draft of my C++ program, I simply used the system() command to run 'samtools view' where needed. I am now wrapping the original program in an R package, and simply running the "system" command does not work. So, I need to run something with the same functionality as "samtools view" from within C++ itself.

@annabelperry I can strongly recommend looking around at the internals of samtools if you need to work with htslib. Samtools is very well put together and I've found scrolling through its code to be a great way to gain insight into how to structure my own programs using htslib.

#include <stdio.h>
#include "htslib/sam.h"
int main(int argc, char *argv[])
{
  samFile *fp = sam_open(argv[1], "r");
  hts_idx_t *idx = sam_index_load(fp, argv[1]);
  bam_hdr_t *h = sam_hdr_read(fp);
  bam1_t *b = bam_init1();
  hts_itr_t *itr = sam_itr_querys(idx, h, argc < 3? "." : argv[2]);
  while (sam_itr_next(fp, itr, b) >= 0) {
    if (b->core.tid < 0) continue;
    printf("%s\t%d\t%d\n", h->target_name[b->core.tid], b->core.pos,
      b->core.pos + bam_cigar2rlen(b->core.n_cigar, bam_get_cigar(b)));
  }
  bam_destroy1(b);
  hts_itr_destroy(itr);
  bam_hdr_destroy(h);
  sam_close(fp);
}

The usage would then be ./a.out alignments.bam chr1:123-456, assuming your BAM file is sorted an indexed. Importantly this is making use of the index via sam_itr_querys() so it's not iterating over the entire file.

So I'm looking for a solution which doesn't require iterating through the whole bam file. Unless I am missing something, this solution requires that you iterate through the whole bam file and add "if-else" statements within the "for" loop to extract out your reads of interest. I am trying to avoid this, as it takes so much time. I'm looking for a solution which searches the binary bam file for reads in specific scaffolds/regions WITHOUT iterating through every read

@annabelperry That isn't reading through the entire file. sam_itr_querys is accepting a string like chr1:123-456 and then making an iterator over just those using the index (idx). This is all that samtools is doing internally. I've added a bit onto the end of this message to hopefully clarify that.

Thank you, @DevonRyan, that clears it up!

Have you considered using a sorted and indexed bam file? e.g.
samtools sort in.bam > out.sorted.bam
samtools index out.sorted.bam

This makes arbitrary region lookups trivial. I am more familiar with pysam (HTSlib for python), which has the AlignmentFile.fetch() functionality.
Surely there is an analog in HTSlib API, though I don't know what to search to get it. It is in samtools as demonstrated here.
This answer seems to put the components together that you are interested in, among them sam_index_load(), sam_itr_querys().

My bam files are already sorted and indexed - I will check out this other answer


text:
Time complexity of the Neighbors function
I implemented the following find_neighbors_with_expected_hamming_distance function that generates all k-mers of Hamming distance at most d from the given Pattern.
from .Hamming_Distance import find_hamming_distance

dp_map_that_stores_all_neighbors_with_expected_hamming_distance = {}

unique_nucleotides = {'A', 'C', 'G', 'T'}

def suffix(pattern):
    if len(pattern) <= 1:
        return None
    else:
        return pattern[1: (len(pattern))]

def first_symbol(pattern):
    if len(pattern) == 0:
        return None
    else:
        return pattern[0]

def find_neighbors_with_expected_hamming_distance(pattern, expected_hamming_distance):
    if expected_hamming_distance == 0:
        return {pattern}

    if len(pattern) == 1:
        return {'A', 'C', 'G', 'T'}

    if (pattern, expected_hamming_distance) in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        return dp_map_that_stores_all_neighbors_with_expected_hamming_distance[(pattern, expected_hamming_distance)]

    neighborhood_set = set()

    suffix_neighbors_set = find_neighbors_with_expected_hamming_distance(suffix(pattern), expected_hamming_distance)

    for suffix_neighbor in suffix_neighbors_set:

        if find_hamming_distance(suffix(pattern), suffix_neighbor) < expected_hamming_distance:

            for nucleotide in unique_nucleotides:
                neighborhood_set.add(nucleotide + suffix_neighbor)

        else:
            neighborhood_set.add(first_symbol(pattern) + suffix_neighbor)

    if (pattern, expected_hamming_distance) not in dp_map_that_stores_all_neighbors_with_expected_hamming_distance.keys():
        dp_map_that_stores_all_neighbors_with_expected_hamming_distance[
            (pattern, expected_hamming_distance)] = neighborhood_set

    return neighborhood_set

I want your suggestion on my time complexity calculation:
For a given pattern of length N and D as maximum expected hamming distance, I calculated time complexity as -
$\sum_{n=1}^N$ $\sum_{d=1}^D$ $n\choose d$ * $3^i$
Reasoning: For any length n, we will choose any d letters (nucleotide in this case) to be replaced with the other 3 nucleotides that is different from the nucleotide present at that particular position.
For example: Pattern = ATGCAT and D = 2, then suppose for one of the case, I select 2nd and 4th positions to be replaced.
Then, at the second position (current nucleotide = 'T'), I can substitute any nucleotide from the set {'A', 'C', 'G'} and
for fourth position (current nucleotide = 'C'), I can substitute any nucleotide from the set {'A', 'T', 'G'}.
Kindly tell me if this approach is correct?
EDIT:
I have provided the pseudocode for the Neighbors function as follows:
Neighbors(Pattern, d)
    if d = 0
        return {Pattern}
    if |Pattern| = 1
        return {A, C, G, T}
    Neighborhood ← an empty set
    SuffixNeighbors ← Neighbors(Suffix(Pattern), d)
    for each string Text from SuffixNeighbors
        if HammingDistance(Suffix(Pattern), Text) < d
            for each nucleotide x
                add x • Text to Neighborhood
        else
            add FirstSymbol(Pattern) • Text to Neighborhood
    return Neighborhood

The symbol • means concatenation of 2 strings. For example, 'S' • 'V' = 'SV'.
If we remove the first symbol of Pattern (denoted FirstSymbol(Pattern)), then we will obtain a
(k − 1)-mer that we denote by Suffix(Pattern).
Illustration:
Given pattern = 'ACG'
Expected Hamming distance (D) = 1
Answer = {'CCG', 'TCG', 'GCG', 'AAG', 'ATG', 'AGG', 'ACA', 'ACC', 'ACT', 'ACG'}
Note that the Answer is a set of all strings who have hamming distance of 1 with respect to the given pattern 'ACG'.

I'm confused by the "expectation" part. Isn't this deterministic?

Also, I'd recommend adding a pseudocode explanation of your approach in addition to the python code.

@Throckmorton, thank you very much for having a look. I have added the pseudocode with an illustration.

Let's first establish a lower bound on any solution:
Given a hamming distance $D$ and an input pattern $P$ of length $N$, there are $\sum_{k=0}^D{N \choose k} 3^k$ other strings within hamming distance $D$ of $P$. Of course in the worst case, when $D=N$, there are $4^n$ possible strings.
We can also get a closed bound on the number of solutions. We can show that the sum of the first $D$ binomial coefficients for a fixed $N$ is bounded by ${N \choose D} {N-(D-1) \over N-(2D-1)}$ as shown here. Therefore we have the following closed upper bound on the number of strings within hamming distance $D$ of some pattern of length $N$
$$
\sum_{k=0}^D{N \choose k} 3^k \leq {N \choose D}  {3^D (N-(D-1)) \over N-(2D-1)} \leq \left(\frac{3N}{D}\right)^D {N-(D-1) \over N-(2D-1)}
$$
However, for your approach, you are going through recursively, so when the subproblem is on strings within $D$ of the suffix of length $n$, you have $\sum_{k=0}^D{n \choose k} 3^k$  of these. Given that you have $N$ subproblems (suffixes), at first glance it does indeed look like you have roughly $\sum_{n=1}^N\sum_{k=0}^D{n \choose k} 3^k$ iterations of the innermost for loop. However, in that for loop, you are also computing the Hamming Distance function, which is $O(n)$. This means that you have roughly
$$
\sum_{n=1}^Nn\sum_{k=0}^D{n \choose k} 3^k
$$
computations.
I would suggest finding a way to get rid of having to compute the hamming function. Since you are constructing the sequences yourself, you can keep track of each sequences hamming distance from the pattern during construction.

Excellent! I missed the complexity for computing Hamming distance. Also, you can correct the first sentence - "Given a hamming distance  D"

Thank you very much @Throckmorton.


text:
Using txt file to query Ensembl API
I am trying to use an external file with a list of gene names to query the ENSEMBL IDs using the perl API. This script works:
#!/usr/local/bin/perl

use strict;
use warnings;
use Bio::Ensembl::Registry;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'useastdb.ensembl.org',
    -port => 5306,
    -user => 'anonymous'
);

my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );

my @genes;
@genes = ('CCL2', 'CCL5');

foreach $name(@genes){
    my $gene = $gene_adaptor -> fetch_by_display_label($name);
    print "GENE ", $gene->stable_id(), ",";
}

However I can't figure out how to get this to work with an external file that contains a list of genes like this:
genes.txt
CCL2
CCL5

Here's what I tried:
#!/usr/local/bin/perl

use strict;
use warnings;
use Bio::Ensembl::Registry;

my $registry = 'Bio::EnsEMBL::Registry';

$registry->load_registry_from_db(
    -host => 'useastdb.ensembl.org',
    -port => 5306,
    -user => 'anonymous'
    );

my $gene_adaptor = $registry->get_adaptor( 'Human', 'Core', 'Gene' );

open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    my $gene = $gene_adaptor -> fetch_by_display_label($_);
    print "GENE ", $gene->stable_id(), ",";
}

But I'm getting this error:
Can't call method "stable_id" on an undefined value at ./get_ensembl_ids.pl line 22, <IDs> line 1.
I was trying to follow this tutorial for reading txt files in perl, but what am I missing for getting perl to actually fetch the ENSEMBL ids with each of the .txt file lines?

Similar to reading files in other languages, it is necessary to trim the newlines from each line. If you add chomp; to the first line of the while loop, the code runs just fine (note that I also print a newline at the end of the loop).
open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    chomp;
    my $gene = $gene_adaptor -> fetch_by_display_label($_);
    print "GENE ", $gene->stable_id(), ",";
}
print '\n';

The new loop produces the following result:
>$ perl query.pl
GENE ENSG00000108691,GENE ENSG00000271503,

Additionally, you may need to catch any exceptions/errors thrown if you have an invalid gene in your input file:
open (IDs, '<', "genes.txt") or die "can't read input file";
while (<IDs>){
    chomp;
    eval {
        my $gene = $gene_adaptor -> fetch_by_display_label($_);
        print "GENE ", $gene->stable_id(), ",";
    };
    if ($@) {
        warn "Error:\t$@\n"
    }
}
print '\n';

Thanks @Throckmorton for the response. Your code is still throwing the error at me, but manages to retrieve the Ensembl ID for CCL5: `Error:  Can't call method "stable_id" on an undefined value at ./get_ensembl_ids.pl line 20,  line 1.
GENE ENSG00000271503,\n% ` Any idea why that might be? Can I ask what version of perl you're using?

I am running 5.26

@user45678789 I will say it is a bit suspicious that you have use Bio::Ensembl::Registry; as well as my $registry = 'Bio::EnsEMBL::Registry';, I would assume the case should be the same in both lines.


text:
Genome annotation
I have helped a lab sequence a mitochondrial genome. I used then MITOS to annotate the genome (warning them that I did it just for curiosity as I am not experienced). They submitted the sequence and annotation to genbank but it was refused for different reasons.
It's the genome of a species for which there is no reference. It's a jumping spider and the genome is almost identical to the one below. I used MITOS web server and I got an annotation almost identical to this one (http://mitos.bioinf.uni-leipzig.de/result.py?hash=bJuyJw4b) which corresponds to this sequence (https://www.ncbi.nlm.nih.gov/nuccore/MH891570 ). But the submitted annotation to Genbank is slightly different so I guess they manually edited it or they used another software. I am not sure I can share more specific information without the permission of the lab.
You can download a GFF file or the raw sequence at the page in the link (for the similar sequence).
My questions are as follows:

after applying MITOS, is there any manual or automatic review that should be applied?
Is it common that the annotation returned by MITOS is not good enough?
Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)
Should we edit the annotation ourselves or there is a software that does it itself?
How often using MITOS is enough or how often manually editing is needed?

The sequence corresponds to an arachnid. The reasons given by Genbank are:

Some or all of the sequences contain incorrect annotation and/or
contain little similarity to other sequences in the database based on
BLAST searches. Please also review your nucleotide sequences for
issues such as low quality data or frameshifts.
Incorrect annotation is:

coding regions lack stop codons
CDS features do not have the correct nucleotide spans or location
for example, nad5 is longer at the amino terminus compared to
similar proteins in the database and extends into the flanking tRNA.

The first part of their reply is apparently generic as I have seen the same answer in another submission in which I didn't participate. I am not sure it's even verified as there are quite good BLAST matches.
I have no idea about the other part. I just applied MITOS.
I used the invertebrate mitochondrial translation table option in MITOS. I guess it's the correct one as the result quite makes sense except for small adjustments.

Any regions annotated as protein sequences must have a complete protein sequence, including start and stop codons. It's somewhat common for mitochondrial genomes to require polyA extension for the stop codons to be completed; I expect that NCBI is asking that the submissions make this explicit. See, for example, our lab's Nippostrongylus brasiliensis mtDNA submission:
https://www.ncbi.nlm.nih.gov/nuccore/NC_033886.1
As an example in this genome annotation, the ND5 gene requires extension, as mentioned in the note:
3858..5439
/gene="ND5"
/locus_tag="B2K98_mgp10"
/note="TAA stop codon is completed by the addition of 3' A
residues to the mRNA"
/codon_start=1
/transl_except=(pos:5439,aa:TERM)
/transl_table=5
/product="NADH dehydrogenase subunit 5"
/protein_id="YP_009347353.1"
/db_xref="GeneID:31081492"

It doesn't surprise me that an automated tool is not perfect. The mitochondrial genome is small enough that manual verification can be feasibly done. You should at least look at the annotated information to make sure it makes sense; the response from NCBI suggests that it doesn't (e.g. nad5 extends into an annotated tRNA sequence).
As a direct / explicit response to your questions, just in case my answers are useful:

after applying MITOS, is there any manual or automatic review that should be applied?

I don't know. You should at least look at the annotated information to make sure it makes sense.

Is it common that the annotation returned by MITOS is not good enough?

I don't know; I have never used MITOS.

Or could that indicate that maybe there is a problem in the assembled genome? (looks unlikely except for a few SNPs as three assemblers return almost the same sequence)

I don't know; you haven't provided the assembled gene sequence or annotation. The response from NCBI suggests that there is a problem with either the sequence or annotation (e.g. nad5 extends into an annotated tRNA sequence).

Should we edit the annotation ourselves or there is a software that does it itself?

If there is an error in the annotation through automated methods, you should fix it manually.

How often using MITOS is enough or how often manually editing is needed?

I don't know; I have never used MITOS.

Thanks for the answer. Unfortunately, it's not clear enough. Let me add what I infer from your answer.

After applying MITOS, is there any manual or automatic review that should be applied?
-> There is no automatic review method and maybe manual review is needed.

Is it common that the annotation returned by MITOS is not good enough? -> You don't have an answer but it would not surprise you if it is.

I am a bit lost with so many edits but I would say the other questions are the same as these two.

You might be luckier by submitting to INSDC by the ENA gate instead of NCBI.
To do so you first have to convert your GFF into the EMBL format used by ENA with the help of EMBLmyGFF3. During the conversion you can use the --tranlate option that can help to see if there is anything wrong in the ORFs (e.g. frameshifts).
Incomplete genes is supposed to be supported by INSDC databases, so I don't get why NCBI complains about it. Either their pipeline does not take care of this, or it is something specific to Mitochondrion annotations. Using EMBLmyGFF3 incomplete CDS will be handled properly and the location field in the resulting EMBL flat file will have the information that it is incomplete (see 3.4.3 Location examples in the documentation http://ftp.ebi.ac.uk/pub/databases/embl/doc//FT_current.html#3.4)
Then once you have your EMBL flat file you can check by yourself if there is anything not accepted by the DB using the ENA webin-cli (https://github.com/enasequence/webin-cli/releases/tag/v3.7.0)
It there are errors they will be listed with details and you will have to fix your file accordingly. When everything is fixed you can submit to ENA. Once accepted by ENA it will be within 24 hours within NCBI.

Thanks, I will comment them that option.


text:
P value recalculation from multiple analysis (Transcriptomics)
I've done a differential expression analysis with transcriptomic information (microarrays and RNAseq) from different datasets by applying for some of them t-test and for the rest Wald-test. I've got many coincidences for these significant genes so I want to recalculate my p-values (or my adjusted p-values, obtained by the BH method) using all these datasets. What is the proper method to do this (if possible)? I am completely new in this field so any additional information would be appreciated.

For cDNA-derived data have a look at the DESeq2 analysis guide.
The general process is to convert all your input data into an integer count matrix (not expression matrix), with genes in rows and samples in columns, and pair it with a metadata data frame that has one line per sample (ordered in the same way as the columns), and indicates the conditions for each sample (including sequencing batch). By incorporating the sequencing batch into the statistical model, it's possible to reduce any technical error associated with consistent linear changes in expression based on the batch (e.g. different runs giving different numbers of total reads).
Microarrays can be combined, bearing in mind that microarray results are often inconsistently presented as normalised / non-normalised, or linear / log. I've found that normalising within-batch then converting to log space usually seems to work well. There is a verbose Users Guide for Limma Voom, but I've found it to be tricky to navigate - it's less friendly to new users.
Combining microarrays and RNASeq is complicated, because the data have different dynamic ranges, and different biases. I think the best thing to do is to process the microarray data separately from the RNASeq data and generate two lists of biologically-interesting genes, then look at the set relationships (e.g. intersection, exclusion).

What is the proper method to do this (if possible)?

Please use dedicated software for RNA-seq and microarray analysis and not naive tests such as the t-test. For microarray one commonly uses limma and for RNA-seq either DESeq2, edgeR or limma-voom, all are available as Bioconductor packages in R.


text:
Search for specific sequence in group of bacterial genomes
I have a group of L. plantarum genomes. I identified the pangenome using Roary. From the pangenome, I was able to identify unique 'genes' for some of the strains. One of these (actually several were, but I'm only interested in one for now) was tagged as 'investigate' in the QC column of the gene_presence_absence matrix, so I did some investigating. Looking at this sequence, I saw it was a hypothetical protein of around 4k bases, and when I ran it through BLAST on the NCBI website, far from being unique, it gave 100% matches to many other L. plantarum strains. Further investigation suggests that this sequence was present on a contig not much larger than the 'gene' its self.
I have two questions.

any ideas why this was identified as unique by Roary but ubiquitous in other members of the species on BLAST (accepting that it's possible it genuinely is unique among my 30 or so strains, but common in the species more generally - it just seems unlikely) - maybe something to do with poorly called open reading frame in prodigal (genome was annotated by prokka)? Or something to do with Roary separating paralogs?

what is the simplest way of confirming if this 4k sequence is/isn't present in my group of 30 genomes?

Any reason you don't just BLASTn your sequence against your 30 genomes to check?

That was my thinking. Would I just use the assemblies to make a blast database and then blast against my 'gene'?

You'd blast your gene against the database, but yes. That will immediately show you if the sequence is indeed unique or not. Then you can focus on the first question.


text:
How to modify dot plot in MUMmer 3 for bacteria comparative genomics?
I was trying to make a dotplot to visualize genome-genome sequence alignment by using the MUMmer 3 software (it is typically used to compare whole genome sequences of bacteria). I runned the same command lines that I wrote in this question post What is the file .mums created by mummer?.
I obtain a plot that is consistent with what I expected but I do not understand why it does not show the axis scales (or labels ? I do not know how to explain but I report below my plot vs the plot made in the MUMmer page example 2.1. mummer so that you can understand what I mean). There could be two solutions:
1. Since my command lines are the same as the one of the example that I found on web, I think that I should modify the mummerplot code to fix this problem.
However, the code is written in C/C++, so not really easy for me. I tried to modify some lines of code but with no success. If you want to have a look at the code and give me some suggestions, it is available on Github.
2. But I wonder also if it is the only way to adjust the mummerplot plot. I found in the MUMmer 3 manual this:

The mummerplot script outputs three files, .gp .fplot .rplot, when run with standard parameters. The first of which is the gnuplot script. This script contains the commands necessary to generate the plot, and refers to the two data files which contain the forward and reverse matches respectively.  [...] Line thickness, color, and many other options can be added or removed from the plot by hand editing the gnuplot script.

Indeed I noticed that my .gp file has some differences with respect to the .gp file of the example that I linked above (the one here). So I could simply modify my .gp file writing it as the one of the example and obtain the plot as I want it, as also the manual suggests. But how can I do it if running mummerplot to generate my plot I obtain the .ps image of the plot and the .gp file at the same time ? The .gp file is not something that I can create before running the mummerplot command...
Below I report, respectively, my plot and the plot of the example (whose aesthetics I would like to imitate):

Thus, I ask if you can give me any hint by answering to one of the two questions or even both or even if there is a third way to solve that I ignore. Thank you in advance !


text:
Pseudocode for gapped BLAST seed extension
I am trying to understand in detail how does the gapped BLAST seed extension DP algorithm works. Despite searching for days now, I could not find anywhere a pseudocode that would allow implementing this in python. Not being a C programmer, it would be impossible for me to grasp this from the BLAST source code.
Can anyone help me with this?
Thank you all. Best,
Daniel

So it turns out the Gapped Blast paper doesn't actually give the extension algorithm. It gives a high level overview and refers to an article in preparation at the time:

The central idea is to consider only cells for which the optimal local alignment score falls no more than $X_g$ below the best alignment score yet found. Starting from a single aligned pair of residues, called the seed, the dynamic programming proceeds both forward and backward through the path graph (Zheng Zhang et al., manuscript in preparation) (Figs 3a and 4). The advantage of this approach is that the region of the path graph explored adapts to the alignment being constructed. The alignment can wander arbitrarily many diagonals away from the seed, but the number of cells expanded on each row tends to remain limited, and may even shrink to zero before a boundary of the path graph is encountered (Fig. 4). The $X_g$ parameter serves a similar function to the band-width parameter of the earlier heuristic, but the region of the path graph it implicitly specifies be explored is in general more productively chosen.

The motivation for this approach is to be a "smarter" version of the banded alignment approach. The issue with banded alignment approaches is that they implicitly won't allow for large insertions/deletions within the alignment, even if the areas around the alignment are high scoring. Increasing the band with to allow for this reduces the speedup gained by the band constraint. This is where the idea of "Alignments Without Low-Scoring Regions" comes into play. This paper is the one cited as "in-prep" by the gapped blast paper, and it describes an early version of the gapped extension approach. However, that  approach seems to have been superseded by the one described here 2 years later by the same lead author (the former paper in '98 having 67 citations and the latter in 2000 having 4700+).
For completeness, I've pasted some pseudocode from the more popular x-drop paper from 2000. Note that $X$ is a user specified parameter, with larger $X$ allowing for a larger search space but slower speed.

Thanks a lot. From what I understand, though, the seed extension is nothing but a "semi-global" alignment, in which the position 0;0 of the matrix is initialized to the score of the kmer match, and the first row and column are initialized to the that score, minus the gap open + gap extend penalties. Right?

Yes, the algorithm is semi-global in the sense that it seeks two prefixes of the input sequences that maximizes the alignment score. However, the initial value of the origin is arbitrary. The algorithm only cares that no cell in the search path falls $X_g$ below the best value seen thus far. The novelty of the X-drop algorithm in Blast is (1) the search termination bounds (i.e. the x-drop) and (2) the greedy algorithm which is faster than typical DP.

I've updated the links to be to pdfs that don't have a paywall. Note that the DP algorithm pictured in the answer is not the main one presented by the paper. The greedy algorithm (which is more efficient) is presented in section 3.

Thanks a lot. It's fun that nowhere it is specified that the seed-and-extend algorithm is a just a case of "semi-global" alignment. You have been superhelpful!


text:
Taiji doesn't run Pagerank
I'm trying to use Taiji to run a pagerank on our combined bulk RNAseq and ATACseq dataset.
Here's how we run it:
taiji run --config config.yml -n 10

Here's config.yml:
input: "input.yml"
output_dir: "output2/"
assembly: "mm10"

Here's a snippet of input.yml (the real one has 8 conditions):
ATAC-Seq:
  - id: cond1_ATAC
    group: cond1
    replicates:
      - rep: 1
        files:
          - path: p1.cond1.narrowPeak
            format: NarrowPeak
      - rep: 2
        files:
          - path: p2.cond1.narrowPeak
            format: NarrowPeak
      - rep: 3
        files:
          - path: p3.cond1.narrowPeak
            format: NarrowPeak

  - id: cond2_ATAC
    group: cond2
    replicates:
      - rep: 1
        files:
          - path: p1.cond2.narrowPeak
            format: NarrowPeak
      - rep: 2
        files:
          - path: p2.cond2.narrowPeak
            format: NarrowPeak
      - rep: 3
        files:
          - path: p3.cond2.narrowPeak
            format: NarrowPeak

RNA-Seq:
  - id: cond1_RNA
    group: cond1
    replicates:
      - rep: 1
        files:
          - path: m1.cond1.tsv
            tags: ['GeneQuant']
      - rep: 2
        files:
          - path: m2.cond1.tsv
            tags: ['GeneQuant']
      - rep: 3
        files:
          - path: m3.cond1.tsv
            tags: ['GeneQuant']

  - id: cond2_RNA
    group: cond2
    replicates:
      - rep: 1
        files:
          - path: m1.cond2.tsv
            tags: ['GeneQuant']
      - rep: 2
        files:
          - path: m2.cond2.tsv
            tags: ['GeneQuant']
      - rep: 3
        files:
          - path: m3.cond2.tsv
            tags: ['GeneQuant']

We've checked that all the files are there. The NarrowPeak files were created by macs2, the tsv files are 2 column tsv's where column 1 is all EntrezID's we were measuring and column 2 is the expression of that gene.
The output of the pipeline is:
ATACseq
|-Download (empty directory
|-TFBS
| |-all narrowpeak files copied over in bed.gz format
| |-9 motif_sites_part.bed7411-?.gz (where ? is [0-8])
|-openChromatin.bed.gz (all peaks in the dataset)
GENOME (contains the mm10 genome)
RNAseq
|-All tsv files copied over
|-expression_profile.tsv (pasted all txv files together)
|-Download (empty directory)

The manual linked here says at the bottom there should also be two other files in the output directory, GeneRanks_PValues.tsv and GeneRanks.tsv. These are the only files we were looking for. Anybody have a clue why we don't see anything? Thanks.


text:
How to run GSEA analysis on R studio using DEG file list generated from scRNAseq analysis
I want to draw the standard plot of GSEA on Rstudio. I have a data frame that consists of a list of DEGs as below.

Actually, I put the DEG list directly to PANTHER that resulted in the GO or pathway analysis data. However, I want to use GSEA.
Meanwhile, I found nice online site and follow their codes for the GSEA analysis(https://bioinformaticsbreakdown.com/how-to-gsea/).
And, here is the my code.
tmp = read.csv("/data/directory/DEGs_HCvsJSS_DC.csv", header = TRUE) %>% filter(Gene.Symbol !="")

GO_file = "/blue/zhangw/my.kim/HCvsJSS/c5.all.v7.4.symbols.gmt" # got this from http://www.gsea-msigdb.org/gsea/msigdb/collections.jsp

gene_list = tmp$avg_log2FC

names(gene_list) = tmp$Gene.Symbol

gene_list = sort(gene_list, decreasing = TRUE)

gene_list = gene_list[!duplicated(names(gene_list))]

head(gene_list)

res = GSEA(gene_list, GO_file, pval = 0.05) # this caused the error as follows. Error in build_Anno(TERM2GENE, TERM2NAME) :  argument "TERM2GENE" is missing, with no default

It's really dissapointing because I found this relatively easy protocol for GSEA, but found the error, lol.
Do you guys have any suggests for this troubles? or any recommendations for GSEA analysis that can be done with DEG list generated from seruat scRNAseq data?

I'd probably just use fgsea. Your code won't have to change much, it also just uses a ranked list.


text:
How to export summary data from Bowtie2 for MultiQC to read?
BowTie2 puts out summary info to the terminal but doesn't allow me to save each to a separate file. How can I go about saving the output of BowTie2 so that MultiQC can read it and create a website from it?

After trying on my own this is the solution I came up with.
My solution was to use a pipe of the stderr into a new file.
bowtie2 --local --threads 4 -x reference_genome.fa -f -U inputdata.fa -S OutputData.sam --no-hd --no-sq --no-unal -k 5 2> 'bowtie2_summary/inputdata.txt'


text:
How to find disease association of disgenet with string db in homo sapiens
I want to work string-db to read the protein-protein interaction from it and I look up disease association from disgenet but the Ids do not seem to match-up and I do not know what to do.
for example for protein INSR which is mostly afflicted with diabetes and similar diseases I found that from 'disgenet' but I do not find the ids in string, is the a convertion table some where?

protein1 protein2 combined_score
9606.ENSP00000000233 9606.ENSP00000272298 490

PS: I think this question may sound stupid to you but I am an absolute noob and would be glad to help me out.


text:
Compressing read data (converting fastq to fastq.gz) on Windows
I downloaded some of my data in fastq format instead of fastq.gz. I want to upload my data on a server now.
I was wondering if there is a way to convert my fastq files to fastq.gz on windows (I read some stuff about python, but I'm not quite sure) before uploading them on the server, since it takes a lot of time uploading the fastq files.

How about downloading onto the server directly? How did you download, as when it is from NCBI it will be faster to just redownload and gzip-compress on the fly.

7zip is file archiver capable of running gzip-style compression, and runs on Windows: https://www.7-zip.org/

If your Windows 10 build number is 17063 or later, you can use tar command to compress file in Windows PowerShell:
tar -vczf output.fq.gz input.fastq

The -z flag means to compress archive with gzip.

Bioinformatics is 99% done in Unix-like environments, and trying to come up with Windows solutions is unnecessary effort. The Windows Subsystem Linux (WSL2) provides you with a native Linux distribution within Windows that easily allow you to do any Unix-ish thing such as running gzip. It is straight-forward to install, see its online documentation.

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.

You could be a little less condescending especially given how late your answer is.


text:
Long-Covid (LC) GWAS project proposal: should I include a non-LC GWAS?
I am working on a project proposal for an exam, it involves a budget of 1M euros and I got to design a genomic study. I chose a GWAS to assess altered gene expression after covid19 that leads to long-COVID. I thought about 3 cohorts:
Case cohort1: Long-Covid patients (N: 3,000)
Case cohort2: Fully recovered patients (N: 3000)
Control cohort: healthy patients. (N: 6000)
Would it make sense to run two GWAS for the two case cohorts against the control cohort (cohort1/controls and cohort2/controls) then do a miami plot to see if there are some hits that differs from the two cohorts, and pick the ones only presented in the Long-covid for further assessment?
The study will also be supported by a validation step for both gwas. What do you think about this idea? Would it make sense or the fully-recovered patients against control cohorts is totally useless? Thank you in advance

In the end I manage to do 2 cohorts
1: Long Covid Patients (1300N)
2: NON Long Covid Patients that had Covid (1300).
It was the most rational thing to do, and professor liked it. So that's how I manage in th end.


text:
How to compare 2 whole genome sequences
I have the files of 2 persons that got whole genome sequencing, one by nebula genomics and the other by dante labs.

I want to compare their genomes to see if they are father and son.
I want to see exactly which genetic sequences they share.
In the shared sequences, I want to know what mutations exist only in the son.

Which program should I use?
Each company provide fastq and vcf files to the genome. so I have each genome in fastq and vcf. which is better to use?
(I have the files on my hard drive, so is it better to use offline program?)

I’m voting to close this question because I feel like this is not a precise technical question and it cannot be answered here.

This seems like a valid question to me; it is similar to bioinformatics projects I've worked on in the past, and I think it's reasonable for people to ask questions regardless of how complex the answers seem to other people.

I expect that answers to this will depend on what the input files look like. Could you please show the first few lines of the fastq files, and the first few lines of the VCF files (including a few lines beyond the header information)?

I think you can't do (2) and (3) without the mother.
(1) is doable to a very high accuracy. It is better to start with fastq. You need to generate a multi-sample VCF first. Here is a high-level overview: a) align reads with bwa-mem; b) sort alignments with samtools and then index the resulting BAM; c) call variants jointly with freebayes (or GATK if you like challenges). After you get the two-sample VCF, you go through every site where the father is homozygous and it is a high-quality site in both samples (e.g. average coverage, high mapping quality, not in difficult regions, etc). Among these sites, check if the child has the paternal allele. If you have a father-son pair, the child is expected to have all the paternal alleles in theory. In practice, if >99% of paternal homozygotes are present in the child, the two samples will almost certainly be a father-son pair. The actual threshold may depend on how you define "high-quality" sites. Errors may lead to apparent mendelian violations.
In principle, it is possible to do a similar analysis with the two provided VCFs, but you need to be very careful about various technical artifacts.

Assuming you've got genotyped VCF files from the same type of data, and they are similarly normalised, I expect that the VCF files will be easier to work with initially.
That's a lot of assumptions, so until there's more information about the input data, I'll just point you to bcftools merge, which I think would be a reasonable first step given those assumptions:
bcftools merge -o merged.vcf.gz A.vcf.gz B.vcf.gz

This will produce a merged VCF file that contains genotypes from both files on each line. To look for similarities and/or differences, you would then inspect each line to check whether the genotypes were the same or different (e.g. same: 1/1 1/1; different: 0/1 1/1).


text:
What is protein secondary structure assignment?
I am absolutely new to Structural Bioinformatics (I only started last week). I am working on the secondary structure assignment/prediction (actually I am not sure) of proteins using ML.
When it comes to protein analysis, most of the books or research papers talk about protein structure Prediction:
https://bip.weizmann.ac.il/course/structbioinfo/pdf/SecondaryStructurePrediction_WIS.pdf
Only a handful of books or articles talk about protein structure Assignments:
https://www.biorxiv.org/content/10.1101/2021.02.02.429329v1
Can anyone tell me what is the basic difference between Assignment and Prediction?
Kindly explain with examples.

Do you have any examples for the use of Protein Structure Assignment? This paper uses the term in the title but they don't particularly distinguish the terminology.

To make it crystal clear (and make than pun):

Assignment: a 3D structure is known and the residues are assigned a secondary structure
Prediction: an algorithm predicts from linear (primary) sequence what the secondary structure may be —often incorrectly. This can be homology based, covariance, deep learning etc. But most commonly just the vector norm of window residue composition, e.g. proline is never in a helix and helices and sheets have different periodicity (two for a sheet and three and a bit for a helix).

What do you mean by 3D structure? Is it a Secondary or Tertiary structure?

Tertiary/quaternary structure —which may be of only part of the protein and with gaps due to technical limitations. But I mean specifically a list of atomic coordinates (eg. PDB file) of that chain, possibly as an oligomer, with a ligands or other chains as opposed to say just the name of the fold, e.g. Rossmann fold.

Answer from @jan-wilamowski, converted from comment:
That second paper seems to distinguish between the known data that is used for ML training (i.e. amino acid sequences are assigned some secondary structure) and the attempt to replicate this natural assignment by an "artificial" prediction.


text:
How to reconstruct a target string from a large set of overlapping reads
I got this question as a homework. Does anyone know how the it can be solved?

Consider a short string $r$ as a read from a string $S$ if $r$ is converted from a substring of $S$ by at most $k$ swaps at random locations, where a swap is an operation that switch two adjacent characters, i.e., string $s_1s_2...s_i s_{i+1}...s_n$ is one swap away from the string $s_1s_2...s_{i+1} s_{i}...s_n$. Two reads are called overlapping if the substrings to which they are converted are overlapping. Given the input of $k << |r|$, devise an algorithm to reconstruct a target string $S$ from a large set of overlapping reads of $S$.

The question doesn't explicitly mention which of the possible superstrings of reads should be returned, so I think it's asking for any superstring.

Can a read be used multiple times?

@Throckmorton You mean like a repeat? I think it cannot be used because there is no unique solution.

Is there supposed to be a unique solution?

@Throckmorton Hmm.. I think there can be multiple solutions as well. I mean given the fact that the number of reads are large it may conclude that there might be a unique solution.

I would venture to guess there isn’t a unique solution, as in genome assembly there usually isn’t

Unless it's not obvious, the problem has no unique solution. An obvious example is a string which uses entirely one character; the problem is then to determine its length, which you can only do statistically.
As you've correctly surmised, this is more or less the short-read de novo assembly problem, with one modification: the "read errors" are all swaps.
The usual way that we solve this problem is:

Pick a length (which is customarily called $k$ in bioinformatics, but the problem confusingly already used that, so we'll use $\kappa$ instead), and extract all $\kappa$-length substrings of the reads. These substrings will, from here on, be referred to as $\kappa$-mers.
Construct a subgraph of the De Bruijn graph, where the $\kappa$-mers are the nodes, and there is an edge between two nodes if two adjacent $\kappa$-mers are connected in any read. Label the edges with the number of times this happens in all reads.
Prune the graph of any edges that have low frequency. In your example, the kinds of errors are known, so you could in fact "correct" the pruned edges, or even correct the reads, based on the higher-consensus path.
An Eulerian superpath through the graph corresponds to the desired string.

I'm not saying that this is the best approach for your homework question, but this is how we do it in practice.
To understand the method in depth, here are some suitable references. Some of these are tool papers, and others explain the data structures and algorithms well.

Pevzner et al, An Eulerian path approach to DNA fragment assembly.
Zerbino and Birney, Velvet: Algorithms for de novo short read assembly using de Bruijn graphs.
Sohn and Nam, The present and future of de novo whole-genome assembly.

Can you please explain what do you mean by "if two adjacent k-mers are connected in any read"?

@Alex The explanation in this paper is quite good (if I do say so myself). https://academic.oup.com/bioinformatics/article/27/4/479/198367?login=true


text:
How to demultiplex MinION fastQ files without Porechop?
I recently got fastQ files from a MinION sequencer. I would like to demultiplex these files.
Is there any other software to demultiplex (no porechop) in or outside Artic protocol? I would like to learn options different from porechop.
Thank you!

A quick google finds https://bioinformatics.stackexchange.com/questions/7139/demultiplex-nanopore-reads-with-custom-barcodes, https://www.protocols.io/view/demultiplexing-nanopore-reads-with-last-xj3fkqn, https://github.com/nanoporetech/qcat

The Guppy basecaller itself now has the ability to internally perform multiplexing. That is currently the preferred option.
If for some reason, you still want a stand-alone application to demultiplex already basecalled reads, then qcat is an option: https://github.com/nanoporetech/qcat

I have written up my own semi-manual method for demultiplexing Nanopore reads. As distinct from the other methods that are commonly used, it searches the entire read for potential barcode sequences, and splits out identified chimeric reads (containing multiple barcode sequences) into a separate file:
https://dx.doi.org/10.17504/protocols.io.7vmhn46


text:
Which GWGAS Tool is 'better', MAGMA or LDSC?
I am currently doing a project regarding gene-based analysis and gene-set analysis given a certain GWAS dataset, I only know 2 well known gene-based analysis tool which is LDSC and MAGMA.
As of now, I'm familiar with how MAGMA works, but not really sure about LDSC, I did read LDSC Wiki briefly, and I am leaning towards MAGMA atm, but it's probably just because I'm more familiar with MAGMA.
So, which one do you think is more robust? or at least what is the pros and cons of each software?

Hi Andreas, do the papers describing the methods talk about the differences? If you already looked into that, you could share with us what do you think is/not different and why one/other should be better for what application. Hope it won't sound too mean, I am just wondering if the question can be more specific.


text:
Adding metadata to Seurat object
I have a Seurat object of 8 patients. I want to add metadata to that so that I have origin of each cell. At the moment UMAP just shows a bunch of cells while I want to color clusters by sample
This is my metadata
> str(head(metadata))
'data.frame':   6 obs. of  3 variables:
 $ barcodes: chr  "AAACAGCCAAACTGCC-5" "AAACAGCCAACAGGAT-1" "AAACAGCCAAGATTCT-1" "AAACAGCCAAGGTATA-1" ...
 $ libcodes: Factor w/ 8 levels "1","2","3","4",..: 5 1 1 1 5 7
 $ samples : chr  "PN0252_0005" "PN0252_0001" "PN0252_0001" "PN0252_0001" ...
>

> pbmc
An object of class Seurat
36601 features across 26473 samples within 1 assay
Active assay: RNA (36601 features, 0 variable features)
>

> dim(metadata)
[1] 26473     3
>

Thanks @angel much appreciated

If you have single-dimension per-cell metadata, and it's arranged identically to the cell order in the Seurat object, I find it easier to use the double bracket notation to add metadata to a Seurat object. For example:
metadata$barcodes -> pbmc[["barcodes"]]
metadata$libcodes -> pbmc[["libcodes"]]
metadata$samples -> pbmc[["samples"]]

Answer from @user438383, converted from comment:
Did you try the AddMetaData function? This post may also be useful.
cluster_letters <- LETTERS[Idents(object = pbmc_small)]
names(cluster_letters) <- colnames(x = pbmc_small)
pbmc_small <- AddMetaData(
  object = pbmc_small,
  metadata = cluster_letters,
  col.name = 'letter.idents'
)


text:
How to show that one docking is better than another with Pymol visually
I have the following docking results:
Dock_1 has better energy score than Dock_2.
What I want to do is to show with Pymol that Dock_1 is better visually.
Here I use hydrogen bond, highlighted with Yellow (Dock_1) and Gray (Dock_2).
But as you can see, with that I can't show that Dock_1 is clearly better.
That's why I'm looking for better alternative than hydrogen bond.

In the last paragraph of this post are some suggested things to check that may be useful: https://bioinformatics.stackexchange.com/questions/16192/how-to-calculate-rmsd-and-binding-energy-from-the-already-docked-ligand-receptor

maybe you can use the "per-residue interaction scores" (wich is available via Schrödinger).
On the other hand,to show if a Docking Pose is really "better" you could show specific interactions rather than the complete big molecule.. e.g. an important salt bridge is missing. (hydrogen bonds are not that important, due to the fact that they can be bridged with water molecules)
Sorry, this is not really a solution, but I am not allowed to comment jet -.-

Totally. To further this, a salt bridge contributes something like -2.5 kcal/mol, while a h-bond -1 kcal/mol. pi-pi bond or S-pi about -1 kcal/mol, but that does not account for the improved solvatation term, which plays a role.  Coordinated waters play a massive part in many protein, e.g. G-protein alpha–beta interaction, ribosome, barnase-barnstar etc., and a water can contribute some -3 kcal/mol.

@MatteoFerla my key concern is to show that visually.  Is there  away I can show salt-bridge, pi-pi bond, S-pi bond with Pymol to show the difference of two docking?

@Lucas thanks. Is there a way to show "per-residue interaction scores" visually in Pymol? Is there away we can systematically identify specific interactions that matters? Any reference on this? I appreciate your reply.

@MatteoFerla Also, if I want to show salt-bridge with Pymol. What's the reasonable residue to choose from each ligand and receptor?

Your E looks like a ∆∆G in kJ/mol, which means that they were docked with a forcefield method. You may have among the files a per-residue split of the forces.
For PyMOL, select salty, ((byres (resn ASP+GLU around 3)) and resn LYS+ARG) or ((byres (resn LYS+ARG around 3)) and resn ASP+GLU) should work for salt bridges btween canonical amino acids.


text:
How to plot a PCOA biplot with OTU loadings as arrows
I'm working on gut microbiota data. I want to perform a biplot showing my 2 enterotypes and the major OTU explaining the variance between my two enterotypes.
Something like :  from https://doi.org/10.1038/srep08397
I'm using the package ade4 to plot the enterotypes which gives me a nice plot but I can't find how to add the OTU loading
My script :
library(ade4)
#dist_jsd_3a = distance
#sample_data(physeq_genus_3a_epi)$cluster_3a_fullname) = enterotypes
#obs.pcoa$li = pcoa result from dudi.pco function

obs.pcoa=dudi.pco(cailliez(dist_jsd_3a), scannf=F, nf=2)
s.class(obs.pcoa$li, fac=as.factor(sample_data(physeq_genus_3a_epi)$cluster_3a_fullname), grid=F, col = c("#0073C2FF","#EFC000FF"))

this code gives me this plot :
I also run this following code using the envfit function from vegan package which allowed me to plot the OTU but the resulting plot is not as pleasant as the one above with the s.class function from ade4 package.
library(vegan)
#dist_jsd_3a = distance
#groups = enterotypes
#pcoa = pcoa result from cmdscale(dis)

dis<-dist_jsd_3a
groups <- sample_data(physeq_genus_3a_epi)$cluster_3a
pcoa <- cmdscale(dis)
efit <- envfit(pcoa, t(otu_table(physeq_genus_3a_epi)))
plot(pcoa, col = c("#0073C2FF", "#EFC000FF")[groups], pch = c(21,24)[groups])
##Ordiellipse##
## I use confidence level of 0.95) ##
ordiellipse(pcoa, groups, col = c("#0073C2FF", "#EFC000FF"),
            kind="sd",lty=c(1), conf = 0.95, alpha = 0.05, lwd = 2)

plot(efit, col = "red", cex = 0.9)

See the result here :
I want to know if anyone already plotted a biplot using ade4 package and have been able to add arrows for OTU. If yes, any tips would be appreciated.


text:
color by clusters and sampled in Seurat
I have a Seurat object
I want to have both cluster numbers and coloured cells by sample names like this figure (from a Nature paper)

I have tried group.by argument in Seurat but only samples come not cluster numbers like below
> DimPlot(pbmc, reduction = "umap",group.by = "samples")
>

[![enter image description here][2]][2]
Any help that to get something like the first figure? Having cluster numbers and coloured by sample names please

Here is a solution that makes use of LabelClusters() from Seurat:
# creating a plot and assigning it to the "plot" variable
# coloring by sample id and turning labels off
plot <-
DimPlot(
  pbmc,
  reduction = "umap",
  group.by = "some_random_sample_id",
  label = FALSE,
)

# here is the trick, adding cluster labels to the "data" variable of the "plot"
# without this, the "plot" object would not have the clustering information but only embeddings and sample id
plot$data$seurat_clusters<- pbmc@meta.data$seurat_clusters

# and here is a sanity check to see if cluster labels are transferred to the "right places"
all(rownames(plot$data) == rownames(pbmc@meta.data))

LabelClusters(plot, id = "seurat_clusters")

Sorry I get this error  > LabelClusters(plot, id = "seurat_clusters")
Error in LabelClusters(plot, id = "seurat_clusters") :
Cannot find variable seurat_clusters in plotting data

The problem is kind of obvious, did you check wheter plot$data contains a column called seurat_clusters?


text:
Getting all variants in a region from ensembl
I want to get all the variants (SNPs) in a given region of the human genome. Is there a way on how to do this with the ensembl-API? It should be a representation similar to this:
There I added dbSNP to the region details. Having this data in a JSON or XML-Format would be really helpful! Here's the link to the image.
Edit: I know that the export-function on ensembl is a thing. But then I only get the rs-Numbers. To get all the info I need I would have to run (potentially thousands) all the rs-Numbers through some additional API. Getting that info in bulk right away would be a huge time-saver.

REST API overlap region endpoint. For example.

Thats exaclty what I need! Thanks Emily!


text:
Can the header line contain arbitrary text in FASTA format?
The following is a valid .fasta file content:
>HSBGPG Human gene for bone gla protein (BGP) GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.
Is this also?
>Arbirary_Name_iJustCameUP_with_and_other_local_identifiers GGCAGATTCCCCCTAGACCCGCCCGCACCATGGTCAGGCATGCCCCTCCTCATCGCTGGGCACAGCCCAGAGGGT ATAAACAGTGCTGGAGGCTGGCGGGGCAGGCCAGCTGAGTCCTGAGCAGCAGCCCAGCGCAGCCACCGAGACACC ATGAG.

NCBI's FASTA format description:
https://www.ncbi.nlm.nih.gov/genbank/fastaformat/
NCBI's BLAST page describing valid FASTA input:
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp
NCBI's SNP page describing FASTA format:
https://www.ncbi.nlm.nih.gov/projects/SNP/snp_legend.cgi?legend=fasta
UniProt's FASTA-header description:
https://www.uniprot.org/help/fasta-headers
Wikipedia reference (indicating format origin):
https://en.wikipedia.org/wiki/FASTA_format
FASTA program (origin of format, per Wikipedia):
https://fasta.bioch.virginia.edu/wrp_fasta/fasta_guide.pdf
Harvard PolyPhen page describing FASTA format:
http://genetics.bwh.harvard.edu/pph/FASTA.html
FASTA file format (from Pacific Biosciences):
https://pacbiofileformats.readthedocs.io/en/3.0/FASTA.html

I am not aware of an official fasta format description. The only constraint that I know of is that there should not be any whitespace after >. Apart from that particular position, whitespaces are allowed in the header line.
So, regarding your particular example, both are valid fasta sequences. When in doubt, you can use SeqIO from Biopython, if you can parse your file with the following code, it  is should be a valid fasta file.
from Bio import SeqIO

with open("example.fasta") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

Edit per @Chris_Rands' comment
The code below does the same as above, meaning that SeqIO.parse() takes care of opening and closing the file.
from Bio import SeqIO

for record in SeqIO.parse("example.fasta", "fasta"):
    print(record.id)

For such a simple format, that's actually quite a few odd edge cases with FASTA (e.g. see https://bioinformatics.stackexchange.com/questions/388/are-there-any-databases-of-templates-for-common-bioinformatic-file-formats). There's nothing official about Biopython's definition of FASTA (and it's subject to change, e.g. https://github.com/biopython/biopython/issues/1814 ) but I agree it's a practical solution

Just a note: anything following a space or tab on a > line is considered comments. Those texts are not part of sequence names.

@user172818 why are you saying that? Yes, some tools parse sequence names that way but not all, that's just a relatively common convention, but as far as I know there is no fasta standard. So what is considered a comment is down to the author of whatever tool you're using to parse the file. Is there some standard I don't know about?

@terdon, here is the link about the "no space after >" constraint: http://genetics.bwh.harvard.edu/pph/FASTA.html . To be honest I did not try if common fasta parsers work with a space or not.

@haci yeah, that's just the format that PolyPhen requires, not a general characteristic of fasta which, since it isn't defined by any standard, doesn't really have any limitations at all. Here's another reputable source giving their own definition: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp. I am racking my brain but I can't remember where I used to see sequences like > foo (with a space). Either it was many years ago or I imagined it, or both. I'm still convinced I remember it but the brain does play tricks so I won't press the point.

@terdon This convention has been widely accepted by almost every mainstream tool for decades. Actually, section 2.2 in the fasta user guide gives an example of fasta: >sequence_name1 and description. I don't know if this example was in the earliest version of fasta, but anyway the original authors have accepted this convention, too. You can download sequences in fasta from NCBI or Ensembl and you will see they clearly put all sorts of auxiliary information after a space. Fastq has a paper which follows the same convention.

@user172818 when you say "fasta user guide", do you mean the original fasta tool? If that already had that convention, then that's a pretty strong argument indeed. As for being an accepted convention, absolutely, no question about that. I just know that I've seen all sorts of stuff in various fasta headers in the wild so I'm wary of saying that anything is standardized in the format.

@terdon Yes, I was referring to this PDF from the fasta software website (section 2.2): https://fasta.bioch.virginia.edu/wrp_fasta/fasta_guide.pdf


text:
Output of allelic association doesn't write the rsID
I'm using plink (1.9b5) to do allelic association. My problem is that my output does not write the rsID (SNP), it just writes a dot.
Output:
CHR SNP  BP   A1  F_A      F_U      A2  P       OR
1   .    1014 A   0.18     0.07     C   0.22    3.03

My script:
plink --bed file.bed --bim file.bim --fam file.fam --make-bed --out output

Can someone please help me see what I'm doing wrong?

Can you please show the first few lines of the .bim file? That will help to establish whether or not there is rs information in there for Plink to use.

Answer from @user438383, converted from comment:
It seems likely your input .bim file probably doesn't have RSIDs either - is that the case?


text:
Retrieval of metadata of a batch of 3953 sequence data from GISAID database
I need to retrieve COVID-19 metadata (such as date of collection of sequence data, location, etc.) for a batch of 3953 GISAID accessions from the GISAID database.
I am not very comfortable with GISAID. I would prefer an NCBI database solution.
Any help would be appreciated.

You can download the full db and and put out anything you want with your favourite scripting language.
To do so, you need to request an account first. Once you log in, scroll down to the bottom of the download options. There are 2 meta-data download options, pick what you need.

You need to request an account first. Once you can login scroll down to the bottom of the download options and there's 2 meta data download options, pick what U need.

Well done, hope U won't mind I added it to your answer.

The NCBI link to the metadata database is here
If I get a moment I'll talk you through it.


text:
GSEA: Computing an Enrichment Score (as described in the original paper)
The following paragraph comes from the GSEA paper that explains how to obtain an enrichment score for a gene set. Note that the enrichment score is computed up to a position i in L. There are 2 things I don’t understand about this:
1- first, the following image illustrating the method seems to suggest that that we score a gene set across the whole list of ranked genes, which conflicts with the idea of “up to a position i” in L
2- second, if we’re only considering genes up to a position i, then how are we going to capture the genes that are underexpressed (appear at the bottom of the ranked list) for phenotype A?

Actually I believe I get it now. i is not a constant. We vary it from i=1 to N and consider the i for which $P_{hit} - P_{miss}$  is largest in magnitude


text:
ensembldb equivalent in python
Is there a python equivalent to ensembldb? I want to get genome coordinates for a transcript (like the function transcriptToGenome) but need to do it in python. I wrote a R script that I'm calling from python but I was hoping there was a more direct way to do this.
Thanks

pyensembl should do the job:
>>> from pyensembl import EnsemblRelease
>>> data = EnsemblRelease(76)
>>> data.transcript_by_id("ENST00000506751")
Transcript(..., contig='5', start=140861224, end=140863521, strand='+', genome='GRCh38')

Thanks so much this is perfect. Do you happen to know if there is anything similar for NCBI?

Maybe https://biopython.org/docs/1.75/api/Bio.Entrez.html?

I don't see a way to map coordinates with Bio.Entrez unfortunately.  Thanks though.


text:
How can I make Snakemake ignore all previously produced/cached output files?
I want to rerun a Snakemake workflow, ignoring all previously produced intermediary outputs.
Snakemake has about 100 different flags. I just want to make it easier to find an answer to this. What's the syntax/CLI flag I need to add?

Over at SO it seemed that people felt these questions were not programmy enough, I hoped that bioinformatics SO may be more tolerant. If you down vote a few questions like this I get an asking ban. Don't care if you close it, the answer is there. But down voting?

I was unaware that you get consequences because of downvotes, I apologize for that. Please make some minor edits to your Q and your A, and I'll be able to remove my downvotes. Sorry again!

If you want Snakemake to ignore the previous output of all rules, add the flag --forceall. Otherwise, if you just want to ignore the output of a specific rule to rerun that rule, add --force:
snakemake --forceall -c1

The command line interface documentation is helpful as a reference: https://snakemake.readthedocs.io/en/stable/executing/cli.html#all-options

Why is this answer downvoted? Does it not answer the question? Could the downvoter please speak up?


text:
Where does Snakemake send output that's generated by shell commands?
I'm wondering how Snakemake deals with output to stdout and stderr that is generated when running a rule with a shell: command.
I ran the following rule and only got informed that there was a non-zero exit code. However, I did not get a stack trace. To get access to the stack trace I had to run the shell command manually:
# Shortened for clarity
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    shell:
        """
        augur parse
        """

Why is that? And can I somehow make Snakemake output more detailed information when something goes wrong to help with debugging?

You usually define a log file like below, and redirect the stdout and/or stderr to it. Stdout is >, stderr is 2> and both of them combined is &>.
rule parse:
    input:
        sequences = rules.join_downloads.output
    output:
        sequences = "data/sequences.fasta"
    log:
        my_output/my_log.log
    shell:
        """
        augur parse 2> {log}
        """


text:
Copying files from multiple directories to a single
I have multiple .GFF files, each in a single directory per annotated whole genome, that I wish to copy into a single directory so that I can perform a core genome alignment.
However the following code does not seem to work:
for pathname in "~/ST_0075/annotations/*/*.gff"; do cp "$pathname" ~/ST_0075/alignment; done

The directory structure is the following:
Home
-ST_0075
---Annotations
-----Sample_001
----------001.gff

The goal is to copy from multiple Sample_XXX subdirectories within Annotations the XXX.gff files into another directory named alignment

Tinkered a bit and realized the solution did not neccistate a for loop.
mkdir ~/ST_0075/alignment; cp ~/ST_0075/annotations/*/*.gff ~/ST_0075/alignment

Exactly, and thanks for posting an answer! Note that if you're using bash, you can make it even simpler if you enable globstar. You can do: shopt -s globstar; cp ~/ST_0075/annotations/**/*.gff ~/ST_0075/alignment and that will recurse an arbitrary number of directories under annotations, so you don't need to count the exact number of */ you need.

Hi @terdon, is it ** that is doing the trick?

Yes, try man bash | grep -A3 'globstar$' to see the details.


text:
10X scRNAseq: Sample mix-up
The student who was working on scRNA seq of KO and WT lines has made a mistake and he mixed both lines and generate the final sequencing data. Now, we are having gene expression data but don't know which is KO and which is WT.
It's a genetic KO. Looking at the presence of absence of the KO gene might not be helpful here as its not deep sequenced.
Cells were extracted from two stable line of animal; one wild type and other Chd2 gene knocked out. As a sequencing output I am having one fastq file in the format 28+10+10+90 . I can generate a count matrix out of it; but I dont know which barcode from WT and with one from KO.
Is there any way to segregate the output based on the KO gene expression?
Thanks

You might be able to detect the genetic KO by looking for CNVs at the relevant loci, e.g. via infercnv (or infercnvpy), or copykat. You might get separate clusters on the UMAP from KO/WT, but it's hard to know

If the samples were from different individual organisms you could try looking for mitochondrial variants that differentiate between them. You should have reasonable coverage of the entire mitochondrial mRNA transcriptome even in 10x datasets.
This tool might help:
https://bioconductor.org/packages/devel/bioc/html/mitoClone2.html

Answer from @haci, converted from comment:
If your mutation of interest is "covered" by the "sample read", you should be able to assign which cell comes from which line. By "sampe read", I mean some 90+ bp read that is used for sequencing the insert (no the read for the cell barcode & UMI).


text:
Converting a matrix array to matrix
I have a matrix in R. When I run class(my.matrix) it returns:
"matrix" "array"

Why do I get "array" when in fact I want a plain matrix? The same class output is returned when running as.matrix(my.matrix)).
Why is that?

class() on a matrix will always return:
"matrix" "array"

Simple example:
> class(as.matrix(1))
[1] "matrix" "array"

This is because an array in R is a vector with multiple dimensions, and the matrix is a special type of array with exactly two dimensions, so every matrix is an array and this is what the class command on a matrix is telling us here, see http://adv-r.had.co.nz/Data-structures.html

From the return value of the class function, it looks like this already is a matrix. It just happens that it is considered by R to be both a matrix and an array.


text:
Trinity assembly from many samples
When you combine samples for de-novo transcriptome assembly with Trinity, do you suggest limiting the number of reads for each sample? I had read one of Matthew MacManes' papers awhile back suggesting 40 million reads was a good number for Trinity assemblies of most metazoans.
However, I have 28 samples across 3 treatments that I want to use for DESeq and so, if I follow the advice at BioStars here, that's 28x40 million reads = 1.1 billion reads going into a Trinity assembler... This seems like a troublesome amount. How is it suggested to combine these samples in situations like this?

Assuming these reads are from the same type of thing - and they should be if you're planning to use DESeq[2] on them - then it's fine to combine all the reads into a single Trinity assembly run.
All recent versions of Trinity will perform in silico normalization by default (based on kmer counts), so you shouldn't need to do a manual subsampling of reads. Read normalisation can be done separately; details on that can be found here:
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Insilico-Normalization
My vague memory is that the kmer counts are only used to get an idea of the total transcriptome size (i.e. for determining the average transcriptome coverage). If you want to get a more complete assembly, then consider only choosing samples that are most likely to have different overall transcript expression (e.g. different tissue types, treatments, etc.). If you've got 28 samples, you shouldn't need to worry about replicates; more variation is better.


text:
A unified database for CNV, SNP, Indel and MSI
I am looking for a database or different databases where I can find information on different gene variants for a gene. As an example if I enter PPARG, I could be able to see SNPs, CNV, InDels and MSI.
In my understanding the information on SNPs can be downloaded from dbSNPs, MSI from ClinVar, CNV & InDels from dbVar. Please correct me if I am wrong.
Is there a database that houses all of this information at under one roof?
Your suggestions are very welcome.
Thank you

To Answer my own question, one can use variation viewer to get all the information :)
DisGeNet is another resource.


text:
How to choose the proper way to analyze the differential binding peaks among different samples for ChIP-seq? DiffBind or CSAW?
I have used the TMM method in CSAW package to normalize the composition bias for my different ChIP-seq samples and also got the normalized bigwig files. Next step I want to analyze the differential binding peaks between several mutant with WT for specific IP (for example, H3K4me1). But I was confused about how to choose the proper way to do the differential binding peaks analysis. What I know now, there exist two basic strategies I could use:
option 1. I use the DiffBind package to do that. The initial input file is bam files and a series of bed files generated by MACS2 call peaks (such as narrowPeak format). In this package, I will generate the normalization factors for the raw bam files before starting the differential binding analysis.
option 2. I could also use the strategy in the CSAW package to filter out the low abundance region and only keep the higher abundance region for multiple tests between different samples. (as shown in the second example mentioned in this paper(https://f1000research.com/articles/4-1080/v1). In this example, I do not need to do peak calling beforehand and I just by calculating the counts by setting the bin width as a narrow or wide range.
Based on these, I have several questions:

since I want to use TMM composition bias normalization to normalize all bam files for specific IP, if I choose the DiffBind package to get the differential binding peaks based on raw bam files, will it cause some argument since in both TMM composition bias normalization in CSAW and DiffBind will generate their own normalization factors.
is the peak calling by MACS2 is an essential step before the differential binding peak analysis? can I omit the peak calling step, just use the generate matrix generated by TMM composition bias normalization to do the differential binding peak analysis?
if choose the DiffBind to analyze differential binding peaks, for the bed files of two biological replicates, should I use the common regions or the pooling regions of two biological replicate to do the analysis?
If I pool the coordinate of peaks of all samples (macs2 callpeak) and calculate the count matrix based on the TMM composition bias-normalized bigwig files, is it possible to import the matrix directly into DiffBind to do the differential binding peak analysis?

I would appreciate it if someone would like to give me some suggestions on that! Many thanks.

Did you get the answers somewhere in the end?
What approaches did you use? Thanks and regards,
Kylie


text:
The confusion of using TPM (transcripts per million)
It is shown that TPM  values are not suitable for DEG analysis but good for within-sample comparison since TPM normalized the gene length. My question is first: if TPM is not suitable for across sample comparison, then why are some heatmaps using TPM values?  In what scenario are TPM values really useful? Should downstream analysis use TPM or raw count matrix?
When I analyzed sc-RNA seq data, the Seurat package uses normalized counts (which is CPM) for DEG, but TPM is not advised for differential testing. I'm not sure what I'm missing in the understanding. It seems very convoluted to me.
I just don't get the point that TPM is commonly used as an input for DEG testing by Seurat (Seurat findmarker function uses "data" slot, which is normalized counts), but TPM is not recommended for DEG testing for bulk.

Use raw counts as much as you can. Add in various relevant factors as covariates in DESeq2. RNAseq metrics have come a long way but are still misused by people because it's more convenient to compare some sort of normalized metric that give out relevant confounding factors, I guess.

Thank you for your comment！Just curious, raw count data does not account for sequencing depth, but my intuitive assumption is that accounting for sequencing depth is necessary for the comparison of gene expression between samples. Why using count data is more appropriate than TPM or other methods? Thank you!

DESeq2 will account for library depth. Comparing anything between two samples is not robust, one should compare between groups of samples using statistically sound methods. You can of course compare TPM but it leaves room for false results.

The naive per-million scaling methods do not properly correct for the compositional bias between samples. This is especially true if the groups you compare are expected to be very different, e.g. different organs, see here demonstrated on some GTEx data: https://www.biostars.org/p/9465851/#9465854
Why do people use flawed metrics? Well, because some are very common. Per-million is easy to compute and does not change when new samples are added. The latter is convenient, and sometimes per-million might be good enough for visualization. I never do it though, I always use normalized (or vst) counts from DESeq2 or edgeR.
For differential analysis of bulk data one commonly uses raw counts which are then normalized internally by the established frameworks such as DESeq2, edgeR or limma-voom.
Single-cell data, if you consider each cell a replicate, often can go with simpler stats as the large sample such allows tests such as the T- or Wilcox test. Again, plain per-million scaling is probably bad, and alternative methods do exist, e.g. the deconvulution/sum-factor strategy in the scran package at Bioconductor. I cannot speak for Seurat, as I do not use it.
For most downstream applications such as visualization and clustering you probably want to use the properly normalized counts on the log2 scale. For something like heatmaps one would scale the data (Z-score) before plotting to show relative differences. I see no application for raw counts beyond feeding them into a DE testing framework that then normalizes them. The only applications that use raw countsbeyond the aforementioned DE testing frameworks that come to my mind are Combat-Seq (batch correction for bulk data) and sctransform (varaince stab/normalization/regression) for single-cell data.

Thank you very much for such a comprehensive comment. Most of the confusion is cleared except one: you mentioned ：” For most downstream applications such as visualization and clustering you probably want to use the properly normalized counts on the log2 scale“, here do you mean log2(rawcount+1) or log2(TPM+1）？which one is more preferred? I have definitely seen both, but for the sake of learning, I am trying to get a deeper understanding. Thanks in advance

For scRNA-seq, you don't want to normalize for gene length because the most popular 10X technology only sequences a the 3' or 5' end of the transcript. Hence CPM, which normalizes for sequencing depth alone, is popular for single cell data.
This is different from your average bulk RNA-seq, where you typically get full length transcripts (unless you're doing something like DRUG-seq).


text:
BEAUti not accepting taxa tip dates
I'm using BEAUti as part of my pipeline and the sequence names include | to separate the metadata fields - with culture date being at the very end.
When I attempt to auto-configure the tip dates on the sequences however, some of the sequences will not update their date(raw value) and if I double-click and manually change the date it will then reverse back to 0.
Any idea what may be preventing me from manually correcting the dates or auto-config to ignore the culture date after the last |?

I found the solution, one of my Seq IDs had a , in one of the meta-data fields, this appears to have thrown off the software and all subsequent Seq IDs could not be auto-configured.
Simply removing the , solved the issue.


text:
What is the tP statistic from ANGSD?
We were looking for genetic diversity in the exon level of an MHC transcript and analysing the data under ANGSD (Analysis of next generation Sequencing Data) written in C++
We used ANGSD to calculate Thetas and Tajima D neutrality tests and one of the thetas was a pairwise theta, however I am unsure what this means and how it differs from Tajima's D?
http://www.popgen.dk/angsd/index.php/Thetas,Tajima,Neutrality_tests
We then divided the tP over length of each exon to find genetic diversity in the exon level. However, as I am very interested in statistics, if someone could explain what the tP statistic means exactly it would be highly appreciated.

While I can't explain it to you directly, you might check out the section on theta estimators in Korneliussen et al. (2013) which is cited on the above ANGSD wiki page. The equation for this estimator in this paper traces back to the original Tajima (1989).

I don't know Theta but Tajima D is a population genetics test that assumes constant sample size. When the number of segregating sites (non-conservative sites in an alignment) significantly differs from pi (Nei's p-distance between a pair of sequence) one possibility is selection given a constant population size.
What Tajima found was cool: under neutrality this relationship conforms to a beta distribution. Just to explain that again so the nucleotide differences between every pair of sequences in an alignment (pi) against the total number of possible mutations conforms to a beta distribution for a constant population size under neutrality.
Its a cool idea - like really cool - but the violations - any level of 'extinction' will skew the distribution, population movement and importantly so could a bias in sampling.
So what you get is a pair-wise matrix at 5% critical probability. If your data entirely conformed to Tajima's D it would be an amazing result, but that is not the result you are looking for. So if Tajima D gave a particular pair <0.05  AND it was supported by other tests THEN (sounds like code) it provides a basis to consider selection as an explanation.
Again I don't know Theta but if both tests were congruent for putative selection I'd say that was reasonable. However, Tajima D alone is not a cool result, its cool theory however but not cool in practice. I hope that explains it?
Note: did any one ever follow this up in trees? Cool if they did because there will be a relationship between the balance of a tree (e.g. Bayesian 'Skyline' plots) and Tajima D.


text:
Remove from Multi-FASTA by Sequence ID
I want to remove sequence VRE32514 – it doesn’t belong and thus is the reason it lacks additional metadata. However I tried implementing this code from a similar question:
awk 'BEGIN{while((getline<">VRE32514")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' *.fasta | grep -v ">VRE32514" *.fasta

With no success and the proposed solution in R using seqinr doesn't seem to work either for myself.
#load package
library(seqinr)
#load file containing sequences
data<-read.fasta("test.fasta")
#create a vector containing species names: these are the last part of the string
vec.names<-unlist(lapply(strsplit(names(data), "|"), function(x)x[length(x)]))
#find names to keep: indices which are not in the species to remove
species.to.remove<-c(">VRE32514")
vec.tokeep<-which(! vec.names %in%  species.to.remove)
#write the final output
write.fasta(sequences=data[vec.tokeep], names=names(data)[vec.tokeep], file.out="output.fasta")

Below is an example of what my Multi-FASTA appears as and the desired output:
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT
>VRE32514
AATCAGATTGTTCGCCATGAGCTGCTGAGAATGATATCATGGAAGGGCGGCATCGAAACA
GGCTTTAAATTAAGTGTAGGCAAGAACTATAAGTTTATTGAAAGGTATATATCCGAGGAT

Desired Output:
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

You can use samtools faidx to first index your file and then retrieve only those sequences you specify. For this grep all the sequence headers you want and then remove the unwanted ones, followed by retrieval of the fasta entries with faidx:
#/ index:
samtools faidx your.fa

#/ run:
grep '^>' test.fa \             # grep names
| grep -v '>VRE32514' \         # remove unwanted
| tr -d '>' \                   # remove ">" from name
| xargs samtools faidx your.fa  # get the entries you want

Output:
>VRE32491|736|PUH-10C|Blood|2016-12-07
ATGAGATCAGAAAAAGAAATGATGGATTTAGTACTTTCTTTAGCAGAACAGGATGAACGT
ATTCGAATTGTGACCCTTGAGGGGTCACGCGCAAATATTAATATACCTAAAGATGAATTT
>VRE32493|1471|PUH-10N|Tissue/Surgical|2016-12-08
CAGGATTATGATATTACATATTTTGTAAGTGATATAGAACCGTTTATATCTAATGATGAC
TGGCTTAATCAATTTGGGAATATAATAATGATGCAAAAGCCGGAGGATATGGAATTATTC
CCACCTGAAGAAAAGGGATTTTCCTATCTTATGCTATTTGATGATTACAATAAAATTGAT
>VRE32503|1471|PUH-11N|Wound|2017-01-05
CTTACCTTATTGCCCTTGGAAGAGTTAGATAATTACCTAAAGGGCGATAAATTAATAAAG
GTTCTAATTGATAAAGATTGTAGAATTAAAAGGGACATAGTTCCGACTGATATAGATTAT
>VRE32504|1471|PUH-EMEP|Blood|2017-01-10
CATGTAAGAAAGCCAAGCGCAAGGGAGTATGATGATTGCTGCAATGAATTTTGGAATGTA
ACACCTTATGTTATTAAAGGATTGTGCCGTAAGGAAATTTTATTTGCTATTGATCATTTT

In R one could use Biostrings:
library(Biostrings)

fasta <- readDNAStringSet("path/to/your.fa")
fasta_new <- fasta[!names(fasta) %in% c("VRE32514")]

#/ save to disk:
writeXStringSet(fasta_new, "path/to/outputfile.fa")

Your solution in R using Biostrings works perfectly, thank you! However I couldn't understand how to get your UNIXsolution to work. When I entergrep '^>' example.fasta | grep -v '>VRE32514' | tr -d '>' | xargs samtools faidx example.fastaI get the following error message: \nnot found in FASTA file, returning empty sequenceood|2016-12-07

VRE32491|736|PUH-10C|Blood|2016-12-07
[faidx] Failed to fetch sequence in VRE32491|736|PUH-10C|Blood|2016-12-07`

@MicroBiostat It works fine on my machine. Maybe something odd wirh your fasta, e.g. empty lines somewhere?

Here's a solution using Python:
from Bio import SeqIO

sequences = SeqIO.parse('test.fasta', 'fasta')
filtered = [seq for seq in sequences if 'VRE32514' not in seq.id]

with open('output.fasta', 'wt') as output:
    SeqIO.write(filtered, output, 'fasta')
```