text:
Variant vs Allele vs SNP
Coming from a CS background. Reading through the wikipedia page, these all sounds like the same thing: Variant, Allele, and SNP.

Variant/Allele/SNP: Some gene locus that differs from the idea human.

For example, if 99% of humans have a T at some locus, and you have a G, C, or A, then you have a Variant/Allele/SNP at that locus.
Is there intuitive way to explain the difference?

Variant = a locus in the genome where there are differences between individuals
Allele = one of the possible bases/sequences that can occur at the variant
SNP = single nucleotide polymorphism, a type of variant. Here the alleles are just single bases, but other variants can be insertions/deletions of one or more bases.

First some basic biology, taking human as an example. Humans have 22 chromosomes (image source):

As you can see (if you squint) in the image, we have two copies of each of the 22. With the exception of the sex chromosomes, X and Y (seen in the bottom right of the image above), you have the same complement of genes on each of the copies of the chromosome. Traditionally, the term "allele" refers to the sequence of one of the two copies.
We inherit one copy of each chromosome from each of our parents. So, if you have inherited a mutation at position $Ψ$ of chromosome 21 from one of your parents, you will be "heterozygous" for that mutation. This means that of the two copies of chromosome 21 in your cells, one will have the mutant allele at position $Ψ$, and the other the wild-type (i.e. "normal", although that's possibly a misnomer) allele.
SNP stands for "single nucleotide polymorphism". This means that it describes a difference of a single nucleotide between the genome you are looking at and the reference genome. For example, if the reference genome has a G at position Ψ of chromosome 21 and you have a T, then this is a SNP.
The term SNP is also commonly used to refer to polymorphic positions. This means specific, single base positions of the genome where variation is common in the human population. Having differences from the reference genome in such polymorphic sites is usually not problematic. The specific use of the term varies a little between subfields. If you're working on population genomics, you are likely to use SNP only for such known polymorphic sites, while others may use SNP to refer to any single nucleotide variant.
Which brings us to variant. This is a blanket term encompassing any and all forms of variation in a genome. SNPs are variants, but so are larger insertions or deletions (indels), copy number variants, large transposed sequences, the list goes on. You can basically think of variant as the catch all term used when you don't want or need to use anything more specific.
To sum up, we have:

Allele: one of the two copies of a genomic sequence present in (diploid) cells.
Variant: a site where the genome you are looking at differs from a reference genome of that species.
SNP: a single nucleotide variant often, but not always, in a position of the genome where a significant proportion of the species' individuals show variation.

Just to add on to previous answers, it is important that allele is one of the possible versions or sequences that can happen at a given position.
Take SNP. Most SNPs have two alleles, meaning that in the population we only see two possible nucleotides at this position. For instance, everybody either has an A or a G at this position. This type of SNPs are called biallelic.
In rare cases, a SNP can have three different alleles, like A/G/C. These are called triallelic. So already you can see that allele is different from SNP. They are actually completely different concepts.
Also, the word allele has been used for a long time. Scientists first used this to refer to a certain version of genetic variation at one place in a genome. This happened even before we know what is DNA. On the contrary, SNP was only discovered after we can sequence human DNA, which is basically in the past 20-30 years. Therefore sometimes it is hard to distinguish the subtle differences between these terms, because they can carry both original meaning based on primitive understand of heredity as well as newer meaning in the era molecular genetics. It might be helpful to think about what these term are referring to in different time and contexts.


text:
Efficiently count reads overlapping a split feature
I have a bam file with reads, and a set of features. The particularity is that some features are "split", being made of several sub-features. For example, here bed1 and bed2 are made of a simple feature, bed3 contains a split feature:
     0        10        20
     |---------|---------|
bam   >>
       >>>>>>>>>>>>
           >>
     |---------|---------|
bed1 ####
bed2             #####
bed3 ####        #####

I want to count the reads overlapping with each (set of) feature, but importantly I want to count only once a read that intersects both features in a set:
feature   count
   bed1       2
   bed2       1
   bed3       2   -> NOT 3

In R, I can do that with a GRangesList where bed1 and bed2 are list elements containing 1 GRange, whereas bed3 is a list element containing 2 GRanges; find the reads that intersect a list element, then use unique() to avoid double-counting reads that intersect two GRanges in the same list element:
suppressPackageStartupMessages(library(GenomicAlignments))
exple_beds <- GRangesList(list(GRanges("I",
                         ranges =  IRanges(start = 0, end = 3),
                         strand = "+",
                         names = "bed1-a"),
                 GRanges("I",
                         ranges = IRanges(start = 12, end = 16),
                         strand = "+",
                         names = "bed2-a"),
                 GRanges("I",
                         ranges = IRanges(start = c(1,12), end = c(3,16)),
                         strand = "+",
                         names = c("bed3-a", "bed3-b")))) |>
  setNames(c("bed1","bed2","bed3"))

exple_bam <- GAlignments(seqnames = Rle(rep("I", 3)),
            strand = Rle(rep(factor("+", levels = c("+","-","*")),3)),
            pos = c(1L, 2L, 6L),
            cigar = c("2M","12M","2M"))

sapply(exple_beds,
       \(.x) findOverlaps(.x, exple_bam) |>
         to() |>
         unique() |>
         length())
#> bed1 bed2 bed3
#>    2    1    2

Created on 2022-10-07 by the reprex package (v2.0.1)
However this gets really slow with a big BAM file.
So, is there

either a more efficient way to do it in R,
or a way to do it with bedtools and friends without creating 100,000 BED files?

For anyone reading this later, I did not find a satisfying answer and ran the above R code on subsampled bams to allow it to run in a reasonable duration (in my case it's "good enough"). An easy improvement is to parallelize the sapply(). For Python aficionados, I expect RSeQC is doing something similar.

source: Stack Exchange

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

sort: votes


text:
In CellRanger output, what are the reads with an xf:i:17 tag?
When running CellRanger on 10x scRNA-Seq data, a BAM file is output with the aligned reads. These reads have an xf tag, described here:

The bits of this tag are interpreted as follows:

1 The read is confidently mapped to the transcriptome
2 This read's barcode, UMI, and feature combination was discarded in favor of a different feature with higher read support
4 This read pair maps to a discordant pair of genes, and is not treated as a UMI count
8 This read is representative for a molecule and is treated as a UMI count
16 This read maps to exactly one feature, and is identical to bit 1 for transcriptomic reads. Notably, this bit is set for a Feature
Barcode read, while bit 1 is not
32 This read was removed by targeted UMI filtering.

In a real dataset, I have this distribution of xf values:
$ samtools view possorted_genome_bam.bam |  egrep -o "xf:i:[0-9]+" | sort | uniq -c
66263599 xf:i:0
99892363 xf:i:17
  245179 xf:i:19
 8283967 xf:i:25

It is clear to me that 0 are unaligned, 19 = 16+2+1 are discarded, 25 = 16+8+1 are kept (and used in the feature_barcode_matrix); but I do not understand what happens to the 17 = 16 +1. They seem to be good reads, they are aligned on genes, yet they are discarded.
One possible explanation from my reading of this paragraph is that the reads with 17 are PCR duplicates, thus ignored, whereas the reads with 19 are in a group matching different genes (several reads with the same BC/UMI but matching different genes). But this doesn't look right when exploring the BAM file manually.

From the documentation at https://kb.10xgenomics.com/hc/en-us/articles/115003646912-How-is-sequencing-saturation-calculated-

All of these samflag 0x400 reads have an xf tag value of 17, which consist of bits 1 and 16. This also means these reads do not have the xf bit of 8, which mark representative reads from a group of duplicates. The converse isn't true though. The xf17 consist mostly of samflag 0x400 duplicate reads but also of samflag nonduplicate reads.

Thank you for the useful answer, however I don't think it's fully answering the question: that suggests that some of these reads, but not all, are duplicates of some sort, but still leaves to me some mystery about what the other ones are, how the duplicates are detected etc

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


text:
recover files after groupTuple
I have a process that emits BAM files with a key containing the sample name:
process source{
  output:
    tuple val(sample), path("sample.bam") into ch

  [...]

}

1st question: I have several files per sample, let's assume 2 for now. I want to merge them together with samtools merge. I can regroup that channel with groupTuple(), but how can I recover them in the next process? Considering they all have the same name:
grouped_ch = ch.groupTuple(by: 0, size: 2)

process merge{
  input:
    grouped_ch

  shell:
  '''
  samtool merge sample.bam sample.bam > !{sample}.bam
  '''
}

Is there perhaps a way in the input: block to say that I expect a list, and name the list components?
2nd question: In practice, I have a variable number of "sample.bam" files for each sample value, does that change the way to define the input:? It also appears I will need to use the built-in groupKey(), is there documentation available? I have read this Github issue, but it would help if a general description was available.

My preference is to avoid working with files with the same name as much as possible. However, sometimes this might be unavoidable. Fortunately, file name collisions can be avoided when working with multiple input files by using the * and ? wildcards:

When a target file name is defined in the input parameter and a collection of files is received by the process, the file name will be appended by a numerical suffix representing its ordinal position in the list.
The target input file name can contain the * and ? wildcards, that can be used to control the name of staged files.

Note that there is also the stageAs path option (see input of type 'path'), which can also take a name pattern (using the above wildcards) to let you use a variable in your script. For example:
process samtools_merge {

    conda 'samtools=1.15.1'
    debug true
    stageInMode 'rellink'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    ls -nl
    """
}

workflow {

    Channel.fromPath( ['./path/to/files/*.bam', './path/to/more_files/*.bam'] ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple(by: 0, size: 2) \
        | samtools_merge
}

Results:
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [elated_montalcini] DSL2 - revision: 1427d64a1d
executor >  local (3)
[83/25654c] process > samtools_merge (3) [100%] 3 of 3 ✔
total 184
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample1.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample1.bam
-rw-r--r-- 1 1000 985 187911 May 25 12:18 sample1.bam

total 128
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample2.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample2.bam
-rw-r--r-- 1 1000 985 129947 May 25 12:18 sample2.bam

total 196
lrwxrwxrwx 1 1000 985     34 May 25 12:18 mysample1.bam -> ../../../path/to/files/sample3.bam
lrwxrwxrwx 1 1000 985     39 May 25 12:18 mysample2.bam -> ../../../path/to/more_files/sample3.bam
-rw-r--r-- 1 1000 985 197225 May 25 12:18 sample3.bam

If each group contains a variable number of files, you can keep the same input declaration (as above) but you will want use a special groupKey object so that the collected values can be streamed as soon as possible. It's not well documented, but here's a simple example that sorts the input BAM files prior to merging. Once all of the 'samtools_sort' jobs are finished for a particular sample, that sample can proceed to merging without us having to wait for all the other 'samtools_sort' jobs to be finished (which would otherwise be the case if the 'size' parameter had not been specified):
process samtools_sort {

    tag { bam.baseName }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam)

    output:
    tuple val(sample), path("${bam.baseName}.sorted.bam")

    """
    samtools sort -o "${bam.baseName}.sorted.bam" "${bam}"
    """
}

process samtools_merge {

    tag { sample }

    conda 'samtools=1.15.1'

    input:
    tuple val(sample), path(bam_files, stageAs: 'mysample*.bam')

    output:
    tuple val(sample), path("${sample}.bam")

    """
    samtools merge "${sample}.bam" ${bam_files}
    """
}

workflow {

    Channel.fromPath( params.bam_files ) \
        | map { bam -> tuple( bam.baseName, bam ) } \
        | groupTuple() \
        | map { sample, files -> tuple( groupKey(sample, files.size()), files ) } \
        | transpose() \
        | set { bam_files }

    samtools_sort( bam_files ) \
        | groupTuple() \
        | samtools_merge \
        | view()
}

Results:
$ nextflow run script.nf -ansi-log false --bam_files './path/to/*/*.bam'
N E X T F L O W  ~  version 22.04.0
Launching `script.nf` [cheesy_swanson] DSL2 - revision: da1a82a059
[04/643129] Submitted process > samtools_sort (sample3)
[1a/93a939] Submitted process > samtools_sort (sample1)
[93/fde0a0] Submitted process > samtools_sort (sample3)
[db/208bac] Submitted process > samtools_sort (sample3)
[8f/f2a231] Submitted process > samtools_sort (sample2)
[f3/8f0f9f] Submitted process > samtools_sort (sample2)
[bb/cca23e] Submitted process > samtools_merge (sample1)
[sample1, /home/steve/testing/work/bb/cca23e016a641b201b7e42f9db74ef/sample1.bam]
[ff/bd3a2a] Submitted process > samtools_merge (sample2)
[sample2, /home/steve/testing/work/ff/bd3a2ac83a0691d2b80b04ef97b755/sample2.bam]
[94/a4f8eb] Submitted process > samtools_merge (sample3)
[sample3, /home/steve/testing/work/94/a4f8eb483602e9335b28f5d40058c8/sample3.bam]

Note that the val specified as input to 'samtools_sort' is an object of class nextflow.extension.GroupKey. In the example above, the input to 'samtools_merge' is also a groupKey - but it doesn't need to be. If you need a string value for whatever reason, you can just call '.toString()' to get it back:
    ...

    samtools_sort( bam_files ) \
        | groupTuple() \
        | map { sample_key, files -> tuple( sample_key.toString(), files ) } \
        | samtools_merge \
        | view()

Thank you so much for this very comprehensive answer! I had totally missed that section of the input documentation. And as ATpoint notes, this is the second time you help me on a Nextflow question, this has been extremely helpful!


text:
scRNA-Seq pseudobulk for marker detection
A standard approach for scRNA-Seq is to partition the single cells into individual clusters, then use a Wilcoxon test to find markers that characterize each cluster (or other statistical methods that consider single cells as replicates).
Recently, there has been a push towards using pseudobulk approaches for "differential state" (DS) analysis, which avoid the pseudoreplication bias of considering single cells as replicates. See e.g. Crowell (2020), Zimmerman (2021), Thurman (2021), Squair (2021), or the OSCA book.
However, all the above papers seem to use pseudobulk methods for DS of a single cluster between two conditions (e.g. drug treatment), or occasionally to compare a specific pair of clusters. I haven't seen any example of using pseudobulk to find markers distinguishing clusters (comparisons of one cluster against all others).
I can think of a few reasons, for example:

these pseudobulk methods can be more computationally intensive,
perhaps the Wilcoxon test is good enough to simply find markers,
pseudobulk methods can only work when each cluster is represented in several samples, which is hard to guarantee for every single cluster in the dataset.

On the other hand it seems to me a pseudobulk approach should yield more relevant markers (genes that are robustly enriched across biological replicates). So, am I missing something? Is there a statistical reason not to use pseudobulk approaches for finding cluster-specific markers?

If you want a statistical reason not to use pseudobulk to do DE for marker finding, I think it is the same reason not to use any DE method for marker finding: the double dipping problem (see Gao and Witten 2020).
Briefly, if you cluster a dataset based on genes (we can safely say that clustering in PCA space is still based on genes), then asking which genes are differentially expressed between clusters is a statistical fallacy.
We could just as easily divide apples in 2 groups according to their shade of red and then ask which RGB channels exhibit a statistically significant difference.
This means that even in the absence of real structure, you will identify very low p-values for your DE test. The OSCA book section you mention links to another section in the advanced book that tells you the same thing with a small practical demonstration generating a 200-dimensional Gaussian and clustering it, then testing for differences in clusters.
This does not mean, however, that markers found this way are meaningless; it simply means that you cannot use a p-value cutoff to define whether a gene is a marker in this scenario, or say that since you have statistically significant markers two clusters are different.
You can still rank these markers by fold change, difference in percentage expression, AUC, entropy, etc as well as p-value. Moreover, markers identified using the traditional methods tend to make biological sense, and statistically rigorous methods are still in development (there was a beautiful talk by Anna Neufeld on an approach involving count splitting, you can read the preprint here).
The test you suggest is interesting, but I don't think you would solve the double dipping problem; this being said it looks more robust and appropriate than using every single cell as a replicate.

Thank you for the useful addition, indeed the pseudobulk approach does not solve the double-dipping problem, and will suffer from the same p-value inflation as a Wilcoxon on the single cells.


text:
BWA doesn't recognize HiC options
I'm trying to run bwa on HiC data. This is the command I'm giving:
bwa mem -5SP -T0 -t4 -R '@RG\tID:\tSM:\tLB:\tPL:ILLUMINA\tPU:none' REFERENCE.fa PAIR1.fq PAIR2.fq

and here's the error message:
mem: illegal option -- 5

Supposedly, the 5SP option is for HiC data (cite), but my BWA isn't recognizing it. Anybody have any ideas? My initial thought is to update BWA to the newest version, but I don't want to mess with our replicatibility since I am new to the lab and don't want to be updating software when I don't need to. The server is a Mac Pro and bwa version is 0.7.12-r1039. Thanks.

You can use conda to install and test the desired version of bwa to test. This will not interfere with whatever version installed on the server.

Your link with those command line options is 1 year old.  Looks like that version of bwa is 6 years old.  You probably should update.


text:
Unexpected error using "stdin" in bedtools intersect as a piped command
From the bedtools intersect man (for version 2.30.0):
-a BAM/BED/GFF/VCF file “A”. Each feature in A is compared to B in search of overlaps.
   Use “stdin” if passing A with a UNIX pipe.

The stdin part, however, does not apply to version 2.27.1. The help on this version gives no advice for piping.
bam mem ref fq1.fastq fq2.fastq | samtools view -b | bedtools intersect -a stdin -b blacklist.bed -v > align_black.sam
Error: unable to open file or unable to determine types for file stdin

- Please ensure that your file is TAB delimited (e.g., cat -t FILE).
- Also ensure that your file has integer chromosome coordinates in the
  expected columns (e.g., cols 2 and 3 for BED).

I also tried:
(map to bam) | bedtools intersect -a "stdin" -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -a -b blacklist.bed -v > align_black.sam
(map to bam) | bedtools intersect -b blacklist.bed -v > align_black.sam

But these don't work either, either for the same or a different reason. Here is a site that tells me to use the above syntax: https://bedtools.readthedocs.io/en/latest/content/example-usage.html
What is the solution to this if I want to pipe everything together for my users?

Here's a way to skip conversion to BAM and pipe between tools, which performs the same set operations:
$ bam mem ref fq1.fastq fq2.fastq | sam2bed - | bedops -e 1 - blacklist.bed > blacklisted_reads.bed

With this toolkit, the hyphen (-) follows a Unix convention of serving as a placeholder for stdin.

While I have chosen ATpoint as the official answer, another possible syntax to use is
bedtools intersect -a stdin -b blacklist.bed -v < `(map to bam)` > align_black.bam

The main issue with this syntax is that it obfuscates the direction of flow for the pipeline, piping step 2 after step 1 is better in the official answer.


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:
Batch Effects that are Close to Collinear with Main Effect
Two years ago, I started working on a project uses both RNAseq and ATACseq. It's supposed to be a simple Healthy Control (HC) vs disease study. The sample collection was done in 2 phases, and it was clear that there was a batch effect that we could adjust for.
However, recently, we received additional metadata. Plotting a PCA plot showed that there was another, larger batch effect that we didn't account for--location of sample donation. There are 4 different locations with the disease samples being from any of the 4, but ALL HC came from only one of the locations.
I reran the count data through DESeq2 with this design formula: phase + location + disease. It didn't fuss about collinearity like it often likes to do and then pooped out a big list of DEG in the RNAseq.
I could probably run the DEG through GSEA to see if the results match the disease's previously known signatures for biological validity, but what statistical worries should I have about this design matrix? What justification would I need for this statistical asymmetry? Thanks.


text:
Targeted hybrid capture bam files showing reads outside target regions in my bed file
Sorry if this is a stupid question I am new to analysing targeted hybrid capture data.
I am analysing my HTS data from a targeted hybrid capture enrichment method and subsequent sequencing on a NextSeq.
Looking at my bam files, I see there are reads that are overlapping into my non target regions.
For example - my .bed file has two target regions in chr17:

However I still have reads overlapping between 7578739-7578918 (following example is position 7578850):

Quite a lot of these reads don't hit any part of the targeted regions in my bed file. Are these mistakes? How have they been pulled down for sequencing when they aren't in the target regions of the baits?

Are you asking why your pull down experiment was not 100.00000% perfect?  You should read the literature or ask your co-workers who have done pulldowns before what % is typical for samples like yours.

It strongly depends on exactly what your protocol is. Note that there are a good number of tools for QCing capture data (here, here).
Off-target reads are a standard observation in capture experiments; recall that it is target enrichment, not target perfect purification, as commenter swbarnes2 suggests.
The real question is whether the proportion of those off-target reads is too high for your application, or otherwise suggests problems. Those depend on details of the technique that you will likely have to figure out yourself.
Here are some tables from the first of those QC papers (supplement), showing on-target proportions for various capture samples (sorry, it would have been a bit of effort to get text of the tables rather than images):

These data suggest that on-target reads are 60%-80% of mapped reads. Note that they also lose a lot of reads in mapping, so we cannot exclude the possibility of mis-mapped garbage reads contributing to off-target mapping.


text:
NCBI blast for exact match of a short sequence
I'm trying to Blast for exact matches to the sequence: 'ATTGNNNNGCAAACCA' in the human transcriptome using NCBI Blast on its 'refseq_rna' database. However, when I do a basic query I get "No significant similarity found."
Note the N's in the middle of the sequence, if I get rid of those N's, yielding the sequence 'AGCGGATTGCAAAGCAAACCA', I get a match in the 3' UTR of the human MeCP2 gene (which is correct). I don't understand why adding the N's here make it not work.
I looked at the help section and I believe this may be due to shorter sequences being less "statistically significant". Statistical significance of output is completely irrelevant to me, so I set 'Expect Value' to 100000000000000. However, I still get the same result.
I also tried to modify the HTML of the page to allow me to submit a word size of 4 to the form. This went through, but I still got "No significant similarity found".
Could someone help me out with this? I feel like searching for a short sequence shouldn't be this hard.

UPDATE: while procrastinating, I worked up the below throwaway code into a more generally applicable (but still bare-bones) script.

ORIGINAL:
I feel like you might want to use something other than BLAST to do this.
BLAST is not really designed to solve this problem, I think it's just easy to use. Also, your N-containing sequence (16nt) is shorter than the version without Ns (21 nt), which can't be correct if you care about "exact" matches.
As a side note, the idea of exact matches is a bit tricky with Ns involved. The whole point of ambiguity codes is that things won't be exact.
And all of this is on top of the short-string problem- most alignment programs won't even return results that don't meet some minimum alignment score, which is exacerbated a bit by the Ns. For how BLAST handles ambiguity codes, you can see section 7.2 here.
Why not use regular expressions instead to search your sequence against the transcriptome?
# python
import re
from Bio import SeqIO

fwd = r'ATTG[ACGT]{4}GCAAACCA'
rev = r'TGGTTTGC[ACGT]{4}CAAT'  # reverse complement

# reference file
seqfile = SeqIO.parse("your_reference_file.fa", "fasta")  # a transcriptome reference fasta

for seq in seqfile:
    # maybe use re.match() to get coords but easier to show it with findall
    a = re.findall(fwd, str(seq.seq))
    b = re.findall(rev, str(seq.seq))
    if len(a) > 0:
        print(seq.id, a)
    if len(b) > 0:
        print(seq.id, b)

When I put this in a file and run it against a test case with ATTGAAAAGCAAACCA, it works:
$ python exact_match.py
the_contig ['ATTGAAAAGCAAACCA']

Alternately, you could just blast the 16 non-N versions of that sequence and combine them, as you've shown.


text:
Finding the relevant annotations for biomarkers
I am looking to find a way to map a list of biomarkers of interest, say genes or bacterial functions, to existing annotation databases. Usually this is done in reverse, as such that based on metagenomics profiling you can see which pathways or genes are present. However, I want to be able to determine which genes or pathways to extract beforehand i.e. as I think it will be quite difficult to condense pathways or genes of interest from annotation data.
Example:
Biomarker gene for better sleep - maps to annotation database at ...
Biomarker gene for glucose metabolism - maps to annotation database at..
.. which in turn allows to directly quantify those biomarkers.
Thanks

Hi @alexander and welcome. How do you propose to 'map', are you looking for a specific software tool/website or else custom code? This question is best described as metabolomics.The question needs more clarity.

Thanks for your comment. For example, I have a metagenomic dataset with annotated genes, based on this, I want to quantify presence of pathways or related genes of interest. However, finding the relevant genes within the annotated file seems to be labour-intensive. As such, I would want to see which of my genes and pathways of interest are present in the annotated file(s). As such, my question is how to 'find' my genes of interest within the annotated files or used database.


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

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

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

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

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

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

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


text:
How does sequence alignment help with protein folding prediction?
I've been seeing using sequence alignment to help with prediction of protein folding, but am unable to understand why. Google returns other research papers for these search terms, hence I'm hoping the stackoverflow community can help thank you.

Where did you see using sequence alignment to help with prediction of protein folding? Perhaps I just too far from the topic to understand the question, but I find hard to understand where the question comes from and I wonder if there is a one correct non-opinion-based answer.

Protein sequences from the same family form a sequence alignment. Some positions in the alignment will be conserved, i.e. the same in all sequences. These are often associated with function, e.g. a catalytic residue at position 50.
However there is a more subtle signal in sequence alignments called covariation. Two positions, say residue 30 and residue 60, may appear to change together across the alignment. We infer that they are close in 3D space in the structure:

Initially statistical techniques (DCA, PSICOV etc.) were used to separate direct from indirect (chained) couplings and extract predicted contacts from sequence alignments. Now deep neural networks have improved performance and even allowed the prediction of residue-residue distances and orientations.
These constraints on the structure allow model building, for example using the CNS software. A good set of constraints effectively guarantees the correct fold. The depth and quality of the sequence alignment is the major factor in the accuracy of the predicted constraints.
This approach to protein structure prediction is the current state of the art. Use of covariation from sequence alignments is the only approach that has reliably given accurate models for larger proteins without templates. See this paper and this paper for the best current methods, and tune into the CASP conference this December to see if anyone can improve on them.

Verified top answer @pohje please mark as "resolved". The technical theory is called "convolution" regarding direct and indirect couplings - you can't just use raw covariancce because you get artefacts, but simplicity you can think of it as just covariance.


text:
Does rRNA depletion protocol give higher number of mapped reads in Intronic regions?
Recently, I have downloaded a publicly available dataset, which are 350 tumor samples. I see the following information from the published paper.

They used Ribo Zero Gold and rRNA was depleted. Strand specific data. After aligning the data I did some alignment quality check with Qualimap RNA-Seq QC tool. I visualised the bam files in IGV. Alignment is good. For all samples 90% alignment rate was seen. I observed that in all samples Higher percentage of mapped reads were originating in Intronic regions. Followed by Exonic and intergenic regions.
I have seen a post here Reads mapped to exonic, intronic and intergenic regions where they say high intronic reads could be because of contamination. I googled about higher reads in intronic regions and found some papers Evaluation of two main RNA-seq approaches for gene quantification in clinical RNA sequencing: polyA+ selection versus rRNA depletion and some other links RIBO-DEPLETION IN RNA-SEQ – WHICH RIBOSOMAL RNA DEPLETION METHOD WORKS BEST? in which they said Greater intronic reads were with rRNA depletion protocol.
And even in this RNA-SEQ tutorial, it is mentioned that - A higher intronic mapping rate is expected for rRNA removal compared to polyA selection.
So, my question:
I am working with lncRNAs. So, I'm using the samples prepared with rRNA depletion protocol. Is this higher intronic rate is common in rRNA depleted dataset or do I have to check anything else to proceed further with these samples?

Hi @maven .. cool username. I'm not really sure this is a bioinformatics question.

It's not so much that you have "intronic contamination" or "genomic contamination", rather you're not selecting explicitly for full-length mature transcripts with rRNA depletion. That is the most common cause for higher intronic read rates. There's nothing you can do about this post-hoc, just continue along.
BTW, many lncRNA's are polyadenylated, so you'll keep them with poly-A selection.

Thanks for the answer Devon. I saw this - Intronic and intergenic mapped reads could be long non coding RNAs, most of them are unannotated [https://www.researchgate.net/post/What_are_the_causes_of_reads_mapped_to_intergenic_region_in_RNA-seq]. My approach is to look for novel lncRNAs.

Yes, it's a non-trivial problem if you can't filter out immature sequences before library prep.

Just to make sure, I also checked the rRNA contamination. And I observed there is very less percentage. 12408836 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 81734 + 0 mapped (0.66% : N/A) 0 + 0 paired in sequencing 0 + 0 read1 0 + 0 read2 0 + 0 properly paired (N/A : N/A) 0 + 0 with itself and mate mapped 0 + 0 singletons (N/A : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)

The read counts I extracted using featureCounts are low compared to total alignments number. Is this because the read counts are from the Exonic regions?

Yes, featureCounts ignore intronic reads by default.

Are there any tools to check DNA contamination?

Probably, you’d want to see how much intergenic expression you have, though you’d need some baseline for comparison.

For all the 350 samples I downloaded...I see 15-30% exonic reads, 70-80% intronic reads an 5-10% interegic reads. so, what % of intergenic tell us there is DNA contamination?

You’d need to mine the literature, but that’s not unreasonable levels given that most of it is probably rRNA and tRNA.

Hi Devon....How to filter our rRNA from the data....And you said more intronic reads could be from immature transcripts....Is it possible to use such data for novel lncRNA detection?

It would be very difficult to use this data for lncRNA detection if you can't be relatively certain that it's high quality. You're not going to have an easy go of it. Regarding rRNA exclusion, blacklist very high expression regions that end up having tRNA or rRNA homology.

I have small basic question...How do I get the sequencing depth of the dataset I downloaded.

Count the number of reads and look at their length distribution. FastQC can aid in this.

I can't possibly see how intron contamination is linked to removal of rRNA depletion.
The only reason it would appear to have increased the number of contaminants is because post-rRNA removal the proporation of intron contaminants has increased against the total remaining RNA content. However, the actual total number of contaminants remains exactly the same pre- and post-rRNA removal. By the same token the RNA of interest lncRNAs will also have proportionally increased, so you get a better depth of its predominance and diversity.
Thats just life and perhaps just filter this bioinformatically.

Hi @Michael thanks for the answer. You mean I should remove these rRNA contamination and then proceed? How do I do that? Any tips please


text:
Problems replicating Kover paper results
I am trying to replicate the results for the creator of this repo:  https://github.com/aldro61/kover2_paper for genotype to phenotype machine learning interpretable (decision tree) prediction.
Specifically in the data section when I try to create the dataset as explained here:
https://github.com/aldro61/kover2_paper/tree/master/data
I make a directory called genomes_dir in the data directory. Then I do:
python download_genomes.py 'mycobacterium tuberculosis/kanamycin/metadata.tsv' genomes_dir

at the end of the download it tells me that it saved file genome_paths.tsv the data directory, not genomes_dir. So I run command:
kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress

and it just hangs with a prompt:
(kover) nobu@gold3forever:~/Desktop/BioInformatics/ResearchAssistant/kover2_paper/data$ kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress
> ^C

it does nothing after hours and hours so I control-C'd to end the process.
Does anyone have any idea what is happening?
I also had to install this to get the datasets from PATRIC:  https://github.com/aldro61/patric_tools

The single quote at 'mycobacterium is not matched so the command prompt is waiting indefinitely for the command to be completed. Maybe you want 'mycobacterium tuberculosis/kanamycin/metadata.tsv' (note closing quote)

Thanks for the catch but now I'm getting some kmer_count = get_kmer_count(tsv_path)   File "/home/nobu/anaconda3/envs/kover/lib/python2.7/site-packages/kover/dataset/create.py", line 136, in get_kmer_count     raise Exception() Exception so I'd better check the data I downloaded.


text:
project ideas for a final year CS undergrad/CS masters student
I am enrolled in an online undergrad CS program in Canada and also Georgia Tech OMSCS where I finished computational photography and AI for Robotics.  I also have some deep learning programming knowledge.
For my capstone project in the undergrad, I was building a nutrition app but I lost interest after finding so many other apps and repos with similar ideas.  So I though since I'm taking an intro bioinformatics course this final term before I graduate, I could do a bioinformatics programming project.
Does anyone have any good ideas?  I'm free labor if I can get approval from my course supervisor for your proposal.  Time length would be approximately 10 weeks.
My languages in order of familiarity: Python, C++, C#, Java, C, Javascript, HTML, CSS.


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

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

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


text:
Seurat Violin Plot: Why do dots align in one row?
I am analyzing a rather small scRNAseq dataset (39797 features across 3522 samples) and I have observed that dots in the violin plot are plotted in one line next to each other rather than in a cloud of dots (see plot below):

I am just wonndering what the reason might be for this?
I have followed the Integration and Label Transfer tutorial using SCTransform (https://satijalab.org/seurat/v3.1/integration.html). After the integration process I switched back to the SCT assay, which should contain the normalized counts. For plotting I used the following command:
VlnPlot(..., assay = "SCT", slot = "data"), again specifying to use the normalized counts.
Can someone explain me why this is the case?
Thank you already in advance!

That is because SCTransform (mostly) outputs discrete values, you can check the @slot of your Seurat object.

Thank you for your answer!Do you know why this is mostly for lower expression values?

It's log-scaled.


text:
How do I obtain a vcf read from a pyvcf.Reader?
I wondered how would I be able to obtain variant based on chrom and location range.
Optimally, this function would be like:
import vcf
vcf_reader = vcf.Reader(filename=some_file_name)
vcf_read = vcf.variants_between(CHROM, START, END)

Does anyone how to do so?

I think vcf_reader.fetch is what you want.  You have to index the file first though.
If starting from regular text VCF you can run:
bgzip example.vcf
tabix example.vcf.gz

Then the Python code is almost the same as yours:
import vcf
vcf_reader = vcf.Reader(filename="example.vcf.gz")
vcf_read = vcf_reader.fetch("NC_023663.1", 7803980, 7804000)
for record in vcf_read:
    print(record)

For this example VCF you get:
Record(CHROM=NC_023663.1, POS=7803981, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803982, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803983, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803984, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803985, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803986, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803987, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803988, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803989, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803990, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803991, REF=A, ALT=[G, <*>])
Record(CHROM=NC_023663.1, POS=7803992, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803993, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803994, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803995, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803996, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803997, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803998, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803999, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7804000, REF=T, ALT=[<*>])

API reference:
https://pyvcf.readthedocs.io/en/latest/API.html#vcf.Reader.fetch
Note the "zero-based, half-open coordinate system" it describes (which is why 7803980 in the second Python argument becomes POS=7803981 from the VCF but then it stops at POS=7804000).


text:
Finding motifs: fasta file with 10,000 sequences
I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that "...the object of type 'FastaIterator' has no len()...". The script however worked when only a single sequence was in the fasta file.
In essence, I need some help figuring out how to:

Prepare an instances object(is this correct?) with microsatellite
motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10
and (AG)10;
Parse my 10,000 long sequence file using this object(?); and
Print an output specifying which of my sequences contain microsatellite motifs.

Appreciate any help I can get on this.
#Seq Utils to find Motifs
>>> pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
>>> sequence = Seq(a)
>>> results = SeqUtils.nt_search(str(a),pattern)
>>> print(results)
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
 56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines

#BrutForce to find Motifs in a fasta file containing 10,000 sequences
>>> def BruteForce(s, t):
...     occurrences = []
...     for i in range(len(s)-len(t)+1):
...         match = True
...         for j in range(len(t)):
...             if s[i+j] != t[j]:
...                 match = False
...                 break
...         if match:
...             occurrences.append(i)
...
...     print(occurrences)
...
...

>>> t = "TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG"
>>> len(t)
34
>>> s = SeqIO.parse("group23.fas", "fasta")

>>> BruteForce(s,t)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "<console>", line 3, in BruteForce
TypeError: object of type 'FastaIterator' has no len()

#BrutForce to find Motifs in a fasta file containing a single sequence
>>> s = SeqIO.read("group23v.fas", "fasta")
>>> BruteForce(s,t)
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]

Bio.SeqIO.parse() returns a SeqRecord iterator. To access the sequence of this sequence record, you need to use the .seq attribute so you should update your s with s.seq.

To expand on haci's answer (and since you mentioned being new to Python) you can loop over the iterator object and use the seq attribute like this:
pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
for record in SeqIO.parse("group23.fas", "fasta"):
    results = SeqUtils.nt_search(str(record.seq), pattern)
    print(results)

(More info on iterators here.)  That will work the same whatever the number of sequences in group23.fas.
I'd be cautious about looking for microsatellites with fixed strings.  If the sequences you're looking for are short and have imperfections among the repeats, you might miss them, depending on how many perfect repeats your search pattern has.  For example with these stub example sequences:
>seq1
TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq2
ACGATGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq3
TGTGTGTGTGTGTGTGTGTGTCTGTGTGTGTGTGTGTGTG

The nt_search output gives you:
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 0]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 4]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG']

That is, one match at position 0 of the first sequence, one at position 4 of the second, and no matches for the third (there's a C I put in the middle).  This might not be relevant for you if you're looking microsatellites much longer than your search patterns, but for shorter ones it might be a problem.

Hi Jesse and @Haci. Thanks so much! I followed both your advice and it worked. I tried something less common, that is, trinucleotide repeats (AGC)10 i.e. pattern=Seq('AGCAGCAGCAGCAGCAGCAGCAGCAGCAGC')
I used the following script

for record in SeqIO.parse("group23.fas", "fasta"):
...     results = SeqUtils.nt_search(str(record.seq), pattern)
...     print(results)
...     print(record.id)

The nt_search output gave the following! (only showing two examples):
NODE_226753_..cov_2.031008 ['AGC...']
NODE_226754_..._cov_2.021041
['AGC...', 265, 268, 271, 274]

Here's a solution that I came up with after help from both haci and Jesse. It involves using the following modules: os, glob, SeqIO, Seq, SeqUtils.
import os, glob

os.chdir("C:\python38\Lib\site-packages\Bio\SeqIO")

from Bio import SeqIO

from Bio.Seq import Seq

from Bio import SeqUtils

patttern = Seq("ACACACACACACACACACAC")

folder_path = 'C:\python38\Lib\site-packages\Bio\SeqIO' #specified the folder where my sequences were stored

fasta_paths = glob.glob(os.path.join(folder_path, '*.fasta')) #glob module finds all the pathnames matching a specified pattern i.e. '*.fasta'

with open("motif4.txt", "w") as f:
... for fasta_path in fasta_paths:
...     print(fasta_paths)
...     for seq_record in SeqIO.parse(fasta_path, "fasta"):
...         results = SeqUtils.nt_search(str(seq_record.seq),pattern)
...         f.write(str(results)) #replace with seq_record.seq to save the list of node names of all sequences that were parsed
...         f.write("\n") #ensures results of the parse through all the sequences in the fasta_path is written to a new line
... f.close()

The nt_search output gave the following output (similar to what Jesse showed):
1 ['ACACACACACACACACACAC']
2 ['ACACACACACACACACACAC']
3 ['ACACACACACACACACACAC']
4 ['ACACACACACACACACACAC']
5 ['ACACACACACACACACACAC']
6 ['ACACACACACACACACACAC', 16636, 16638, 16640, 16642, 16644, 16646, 16648, 16650, 16652, 16654, 16656, 16658, 16660, 16662, 16664]
7 ['ACACACACACACACACACAC']
8 ['ACACACACACACACACACAC']
9 ['ACACACACACACACACACAC']
10 ['ACACACACACACACACACAC']
11 ['ACACACACACACACACACAC', 9174]

It would miss sequences that did not match (AC)10 or those with imperfect repeat motifs. However, in cases where the motif was > (AC)10, nt_search listed down all locations where the motif was found sequentially. For example, in line 6, all locations were listed and showed the motif was in fact 36 nucleotides long in total (pos #16,636 to #16,673).
You can then import the txt. files into a databasing app and analyse the results.


text:
Cobratoolbox unable to identify gurobi solver when passing initCobraToolbox
I am running into this issue lately when I am trying to install the Cobratoolbox on matlab along with having installed the prerequisite gurobi solver and adding it to the MATLAB path. I have run the gurobi_setup.m file as well. However I am unable to uninstall the cobratoolbox, as everytime i run initCobraToolbox.m it fails as it not able to recognize the solver. Is there anything I am missing in the installation? I have attached a screenshot for reference

In my opinion it is better to redirect the question to cobra-toolbox google group https://groups.google.com/g/cobra-toolbox

Not that I would know this package, but Installation should be clear in the README file and if you run into troubles while following it, it's good enough to open an issue on the platform where they develop the package. Not that I would think this question does not belong here, but I think you have much better chances there.


text:
How can I obtain a list of all NCBI gene ID's along with their full name, symbol, and also known as symbols?
I would like to do two things:
(1) Obtain a list of all NCBI (Entrez) gene ID's.
(2) Obtain the "Official Symbol", "Official Full Name", and "Also known as" fields from the NCBI gene pages for each of those genes.
This information, since it is displayed on the NCBI website, must be located in some databases, but I have not had luck finding a tool that will allow me to extract this information.

All of these data are in the gene_info.gz file on the NCBI FTP site: https://ftp.ncbi.nlm.nih.gov/gene/DATA/
Specifically, you may want to look at the gene_info.gz file. It has the following fields (shown for an example gene below):
                              #tax_id [  1]: 9606
                               GeneID [  2]: 5768
                               Symbol [  3]: QSOX1
                             LocusTag [  4]: -
                             Synonyms [  5]: Q6|QSCN6
                              dbXrefs [  6]: MIM:603120|HGNC:HGNC:9756|Ensembl:ENSG00000116260
                           chromosome [  7]: 1
                         map_location [  8]: 1q25.2
                          description [  9]: quiescin sulfhydryl oxidase 1
                         type_of_gene [ 10]: protein-coding
   Symbol_from_nomenclature_authority [ 11]: QSOX1
Full_name_from_nomenclature_authority [ 12]: quiescin sulfhydryl oxidase 1
                  Nomenclature_status [ 13]: O
                   Other_designations [ 14]: sulfhydryl oxidase 1|quiescin Q6 sulfhydryl oxidase 1|testis tissue sperm-binding protein Li 62n|thiol oxidase 1
                    Modification_date [ 15]: 20200604
                         Feature_type [ 16]: -

If you are starting  with a set of NCBI GeneIDs and want this sort of information for just those, you can use Entrez Direct for that. If you are starting with RefSeq transcript or protein accessions, be sure to check out the other files in the FTP path mentioned above.

I had tried using this file, but a bug in my code made me think it didn't have what I wanted.  Your answer made me look into it again, and I found the bug.  It's exactly what I needed.

Go to Ensembl (https://www.ensembl.org/index.html), click BioMart on the menu at the top of the page. Select the dataset you want to use, then go to Attributes, select the feature of each gene you want to get. Then click "Results" button to make sure that the table looks like what you want. Then output the results to a file.

This is almost exactly what I'm looking for.  The only thing it is missing is as comprehensive a "Also known as" list.  For example, NCBI lists the PPP4C gene has having a total of 7 symbols, while this data set has 2.


text:
Can SBOL encode DNA Assembly linkers?
Where in an SBOL file could you specify the linkers that would be placed as prefixes / suffixes between parts for assembly? Lets say you have a final construct of your design and want to assemble the component parts - could you encode what linkers you used between each component, ie as an annotation? Would this be different between SBOL versions?

Every version of SBOL allows you to specify prefixes/suffixes as annotations on a sequence.

In SBOL 3, the construct would be a Component of type DNA (https://identifiers.org/SBO:0000251), and the linker a SequenceFeature on that Component. The role of the sequence feature might depend on the particular assembly you are describing, and whether you're talking before or after assembly (in those cases where the process changes the sequence), but a good one in many cases would be Nucleotide Binding Site (https://identifiers.org/SO:0001655)
In SBOL 2, the representation is the same, except the construct would be a ComponentDefinition and the linker a SequenceAnnotation.
In SBOL 1, the representation is still the same, except the construct would be a DnaComponent and the linker an SBOL 1 SequenceAnnotation (which has different fields than in SBOL 2).

If the assembly method does not change the sequence, then all three can represent this pretty much equally well.  If you want to represent a sequence that is modified by assembly or the plasmid carriers from which parts are extracted, then SBOL 3 is the best because it cleanly represents deletions during assembly via the sourceLocation field.


text:
Seurat analysis only with matrix with no barcode and features
I have downloaded a dataset from GEO and only the matrix (mx) files are available. How can I perform the analysis without the barcode and the gene features?

what do you mean by gene features? Does the matrix have genes as rownames? Can you share the GEO number? you need gene names or metadata ...

Seurat has different functions that can read different types of data. Please share a few rows and columns of your data and we can try to point you to the right direction (function).


text:
Plot information stored in dataframe within a tree (ggtree)
Hello to the entire Stackoverflow community!
I'm writing to you because I'm currently building a phylogeny with ggtree and I have in parallel a table like this one:
df :
Node_number Number
50   31
1    15
13   25
35   12
40   4
5    2

and a tree such as :

and the idea is that I woule like to incorporate the df$Number information in each Node in the tree and get somethhing like I put in the picture (expected part)
tree data :
structure(list(edge = structure(c(31L, 32L, 33L, 34L, 34L, 35L,
35L, 33L, 32L, 36L, 37L, 37L, 36L, 38L, 39L, 40L, 40L, 39L, 41L,
42L, 42L, 43L, 43L, 41L, 44L, 44L, 45L, 45L, 46L, 46L, 47L, 48L,
48L, 47L, 38L, 49L, 50L, 51L, 51L, 50L, 52L, 53L, 53L, 54L, 54L,
52L, 49L, 55L, 55L, 31L, 56L, 57L, 57L, 58L, 58L, 56L, 59L, 59L,
32L, 33L, 34L, 1L, 35L, 2L, 3L, 4L, 36L, 37L, 5L, 6L, 38L, 39L,
40L, 7L, 8L, 41L, 42L, 9L, 43L, 10L, 11L, 44L, 12L, 45L, 13L,
46L, 14L, 47L, 48L, 15L, 16L, 17L, 49L, 50L, 51L, 18L, 19L, 52L,
53L, 20L, 54L, 21L, 22L, 23L, 55L, 24L, 25L, 56L, 57L, 26L, 58L,
27L, 28L, 59L, 29L, 30L), .Dim = c(58L, 2L)), tip.label = c("t5",
"t9", "t25", "t14", "t30", "t26", "t2", "t10", "t12", "t15",
"t1", "t3", "t6", "t23", "t17", "t18", "t4", "t22", "t27", "t16",
"t29", "t7", "t21", "t20", "t13", "t11", "t19", "t24", "t8",
"t28"), edge.length = c(0.912875924259424, 0.293603372760117,
0.459065726259723, 0.332394674187526, 0.65087046707049, 0.258016780717298,
0.478545248275623, 0.766310670645908, 0.0842469143681228, 0.875321330036968,
0.339072937844321, 0.839440350187942, 0.34668348915875, 0.333774930797517,
0.476351245073602, 0.892198335845023, 0.864339470630512, 0.389989543473348,
0.777320698834956, 0.960617997217923, 0.434659484773874, 0.712514678714797,
0.399994368897751, 0.325352151878178, 0.757087148027495, 0.202692255144939,
0.711121222469956, 0.121691921027377, 0.245488513959572, 0.14330437942408,
0.239629415096715, 0.0589343772735447, 0.642288258532062, 0.876269212691113,
0.778914677444845, 0.79730882588774, 0.455274453619495, 0.410084082046524,
0.810870242770761, 0.604933290276676, 0.654723928077146, 0.353197271935642,
0.270260145887733, 0.99268406117335, 0.633493264438584, 0.213208135217428,
0.129372348077595, 0.478118034312502, 0.924074469832703, 0.59876096714288,
0.976170694921166, 0.731792511884123, 0.356726912083104, 0.431473690550774,
0.148211560677737, 0.0130775754805654, 0.715566066093743, 0.103184235747904
), Nnode = 29L), class = "phylo", order = "cladewise")

does someone have an idea ?

I will have a look @chippycentra later in the week. The idea is cool. I don't really do R, but I will look into it. There are pop gen solutions to this issue but don't think that quite fits here.

You can accomplish this using the ggtree package available on Bioconductor.
First you will need to combine your tree with the data.
library(tidyverse)
library(ggtree)

ftree <- tree$edge %>%
   as_tibble() %>%
   mutate(Node_number = 1:n()) %>% # finds edge numbering
   right_join(data, by = "Node_number") %>% # find internal node associated with your node numbering
   left_join(fortify(tree, ladderize = FALSE), ., by = c("node" = "V2")) # merge tree and data

Note that this will produce a warning, but it does work.
Then plot the tree using ggtree and its geom_nodepoint layer.
ggtree(ftree, ladderize = FALSE) + # plots tree
   geom_tiplab() + # adds tip labels
   geom_nodepoint(aes(subset = !is.na(Number), size = Number, x = branch)) + # adds your data as points
   scale_size_continuous(range = c(4, 12)) # scales point sizes

You can adjust the size of the points with the scale_size_continuous function.
Data
data <- read_csv('Node_number,Number
50,31
1,15
13,25
35,12
40,4
5,2')

and tree is the same as you provided.
Ggtree also has a decent manual available at https://yulab-smu.top/treedata-book/.

thank you very much for your help :-)


text:
How to find closely related genes for a specific gene from WGCNA modules?
I have used WGCNA by integrating several datasets which are processed in a similar way. Identified 17 modules, followed by performed pathway analysis with the genes present in those modules.
Among all the protein-coding genes,ncRNAs and other genes, I'm particularly interested in a specific lncRNA MALAT1.
I found that MALAT1 is co-expressed with more than 800 genes in a module named blue. I have few questions, can anyone please help me in this?
Questions:

I would like to check which are the closest co-expressed genes to MALAT1. How can I know that?

Once I get the closest genes to MALAT1, I would like to make a network with those genes and MALAT1. How to do this network?

The simplest way is to re-calculate (or load, if saved) the TOM matrix and find the genes with the highest TOM to MALAT1. If you have enough RAM to carry out WGCNA in a single block, you can use TOMsimilarityFromExpr, make sure you give it the same arguments as the ones you used in network construction, and select the appropriate column from the result. Then order the column and retain top genes, as many as you think is reasonable. You can also use adjacency instead of TOM, though TOM is a more natural quantity to use if you used in your network construction. Once you figure out which genes you want to retain, you can subset the TOM matrix (both rows and columns) to those genes and you have the sub-network among these genes.

thanks a lot Peter

small question again...I'm interested in modules blue and brown...So, extracted the expression for the genes in those modules and used this way TOM2 = TOMsimilarityFromExpr(t(expr_of_interest),power = softPower,networkType = "unsigned", TOMType = "unsigned") And this gave me correlation values. I can see all correlation values are above 0 only....there are no negative correlation values....why is it like that?

As the function name says, it returns a topological overlap matrix (TOM) which is by definition non-negative. If you want correlation, use cor.

thanks peter. how to get the correlation values for the genes in specific interesting modules or genes from a couple of modules?

Well, you can select the genes using labels, something like (using the notation from WGCNA tutorials) select = moduleColors %in% c(module1, module2,...) (replace module1, module2, .... with the module labels you want to select), then use cor(datExpr[, select]).


text:
classifying samples by TCGA signature
I have some RNA-seq samples from multiple glioblastoma tumours that I'm now trying to classify according to a specific gene signature (from Verhaak et al., 2010) using R. The gene signature is reported as a gene list with specific centroids for each of the 4 clusters (https://api.gdc.cancer.gov/data/941f81a1-05d7-4f84-80ec-534b8dc1ebac). I'm wondering how I can use this signature to classify my samples in R? Would it involve some sort of k-nearest neighbours method?
Additionally, the signature was identified using microarray data, but I am classifying RNA-seq data. Is there any sort of adjustment I should make to the signature to account for this?
Thanks in advance!

Hi inteeresting post and right site to post it :-)  Coudl you summarise the genetic signature here rather than use a weblink please?

Hi Michael! Thanks for the response. The authors of the paper analyzed 200 glioblastoma samples, and used consensus average linkage hierarchical clustering to identify 4 clusters in their dataset. Next, they used ClaNC, a nearest centroid-based classifier that balances the number of genes per class, and identified signature genes for all four subtypes. The gene signature is a collection of 840 genes (210 genes per subtype). I'm not familiar with gene signatures, but for each gene they have 4 centroid values corresponding to each of the subtypes.

This is not a very simple problem. I saved the table and "test.txt" and visualize it below:
library(pheatmap)
x = read.delim("test.txt")

head(x)
  Gene.Symbol            Composite_chr_coords  Proneural      Neural
1      CDKN1B   [36.1:12:12761576-12766570:+] 0.03353670  0.03353670
2       EPB41    [36.1:1:29086215-29319545:+] 0.03882308  0.03882308
3        CLGN  [36.1:4:141529057-141568232:-] 2.12796056  0.41999847
4      PDE10A  [36.1:6:165664525-165995575:-] 1.66012916  0.25349464
5     RALGPS2  [36.1:1:176960923-177152465:+] 0.78971500 -0.01945154
6        TAF5 [36.1:10:105117714-105138812:+] 1.08354870  0.09461930
    Classical Mesenchymal
1  0.03353670 -0.84379400
2 -0.84821990  0.03882308
3  0.41999847  0.41999847
4  0.25349464  0.25349464
5 -0.01945154 -0.01945154
6  0.09461930  0.09461930

You can see there are no so called markers genes and if you look at the scale of the data, it was normalized in microarray and it will be very hard for you to get your RNA-seq to fall into that range, and apply the medoids to identify your clusters or type.
Since you did not provide any data, and it is not trivial to simulate data for this, I can suggest the approach in words here:

Subset your count table to these exact set of 840 genes.

Normalize the data and do a PCA or hierarchical clustering to see if you indeed observe 4 or 3 or 2 subtypes

If you indeed see 4 or 3 clusters this would be straight forward, take the centres of means of these clusters, and correlate with each of the sub types here. Chances are the subtype that is most correlated with one cluster, is the subtype you are interested in.


text:
How to use EPIC DNAm files (.idat) to infer copy number alteration
I have access to EPIC DNA methylation profiles from approximately 40 samples. I want to use these .idat files in order to predict copy number alterations within specific promoter regions of interest. Can you recommend any R packages that could perform this task?
Through research I have stumbled upon epicopy (https://github.com/sean-cho/Epicopy/blob/master/vignettes/Epicopy.Rmd), but I'm uncertain whether this is the best package for this problem.
Cheers,
GB

This is a good question and very interesting area of investigation. The package you sourced looks okay from my 'outside the field' perspective, but I understand your concerns.

I am not too familiar with this, but conumee (https://bioconductor.org/packages/release/bioc/vignettes/conumee/inst/doc/conumee.html) seems promising


text:
How do I Read a Sam File into a String Using ifstream in C++?
I am trying to read a .sam file into a single string in C++. My end goal is to read this string into a vector, where tabs indicate separations between elements of the vector. After making a vector containing all of the data from the sam file, I will separate the relevant data types (chromosome name, locus, etc.) into their own vectors based on position in the master vector.
The sam file was pre-fetched from the SRA database, split into paired-end reads, and aligned to the reference genome using the following Linux commands:
prefetch --output-directory ./ SRR
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR.sra
bwa mem -M reference_genome.fasta SRR_1.fastq SRR_2.fastq > SRR.sam

To make the file easier to read in C++, I altered the .sam file in Linux by removing the header, replacing newline characters with tabs, and replacing gaps between "columns" with true tabs rather than spaces. Since the number of spaces between each column was variable, I accomplished this last part by replacing successively smaller numbers of spaces (12, 11, 10, etc) with tabs. Here is the Linux command line code for the three steps listed above:
grep -v '^@' SRR.sam > headerless_SRR.sam

# Condense sam files into a single line
sed 's/\n/\t/g' headerless_SRR.sam > condensed_SRR.sam

# Replace all spaces in sam files with tabs. Since some columns have multiple spaces, the spaces must be replaced in progressively smaller chunks
sed 's/            /\t/g' condensed_SRR791885.sam > tch_SRR791885.sam     # "tcd" is tabbed, condensed, headerless
sed 's/           /\t/g' tch_SRR.sam > tch2_SRR791885.sam     # Repeat this line with progressively smaller numbers of spaces until all spaces are replaced with single tabs

To open the resulting file and read it into a string, I created the following C++ script, compiled it using the g++ command in Linux, and ran the resulting .out file using "./"
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>
#include <string>

int main(){
    // Creates an input file stream opbject and fills it with data from the sam file
    std::ifstream sam_file;
    sam_file.open("Sample_Sam.txt", std::ifstream::in);
    // Reads the first character in the sam file and stores it as a "character" object
    char sam_character = sam_file.get();
    // Creates a string into which characters can be pushed
    std::string sam_string;
    // While characters are still detected in the sam file, add the characters to a string and move to the next character in the file
    while(sam_file.good()){
        sam_string.push_back(sam_character);
        sam_character = sam_file.get();
    }
    sam_file.close();
}

When I run the C++ script on a dummy sam file containing random blocks of data, it runs fine and I can even output the resulting string to "std::cout" or another file. However, when I run it on a "tch_SRR" sam file, the "./" file hangs without creating an output or anything. Does anyone know how I could fix this?

Are you sure it hangs, or are you instead noticing how incredibly inefficient it is (at least read a line at a time and don't modify the file. SAM files are formatted how they are for a reason)? 2. This seems more appropriate for the codereview stackexchange.


text:
How do I include repeat purity, default slippage, default stutter, and minimum flanking (left and right) in Tandem Repeat Finder's output?
I am attempting to create a markerInfoFile for use in a program called popSTR (GitHub Documentation: https://github.com/DecodeGenetics/popSTR). The marker info file should contain information about microsatellites in a reference genome. Microsatellites are stretches of DNA where a short motif (i.e. "AGG") is repeated multiple times (i.e. three repeats of "AGG" is "AGGAGGAGG"). Here is an example line from the popSTR author's marker info file:
chr10 10589 10601 GCCC 3.2 NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCACCCTTCTAACTGGACTCTGACCCTGATTGTTGAGGGCTGCAAAGAGGAAGAATTTTATTTACCGTCGCTGTGGCCCCGAGTTGTCCCAAAGCGAGGTAATGCCCGCAAGGTCTGTGCTGATCAGGACGCAGCTCTGCCTTCGGGGTGCCCCTGGACT GGTCTGTGCTGAGGAGAACGCTGCTCCGCCTCCGCGGTACTCCGGACATATGTGCAGAGAAGAACGCAGCTGCGCCCTCGCCATGCTCTGCGAGTCTCTGCTGATGAGAACACAGCTTCACTTTCGCAAAGGCGCAGCGCCGGCGCAGGCGCGGAGGGGCGCGCAGCGCCGGCGCAGGCGCGGAGGGGCGCGCCCGAACCCGAACCCTAATGCCGTCATAAGAGCCCTAGGGAGACCTTAGGGAACAAGCATTAAACTGACACTCGATTCTGTAGCCGGCTCTGCCAAGAGACATGGCGTTGCGGTGATATGAGGGCAGGGGTCATGGAAGAAAGCCTTCTGGTTTTAGACCCACAGGAAGATCTGTGACGCGCTCTTGGGTAGAGCACACGTTGCTGGGCGTGCGCTTGAAAAGAGCCTAAGAAGAGGGGGCGTCTGGAAGGAACCGCAACGCCAAGGGAGGGTGTCCAGCCTTCCCGCTTCAACACCTGGACACATTCTGGAAAGTTTCCTAAGAAAGCCAGAAAAATAATTTAAAAAAAAATCCAGAGGCCAGACGGGCTAATGGGGCTTTACTGCGACTATCTGGCTTAATCCTCCAAACAACCTTGCCATACCAGCCCATCAGTCCTCTGAGACAGGTGAAGAACCTGAGGTCGCAGGAGGACACCCAGAAGGTCCAGAGAGAGCCTCCTAGGCCCCCCACCTCCCCCCGTGGCAGCTCCAACCCCAGCTTTTTCACTAGTAAGGCAGTCGGGCCCCTGGGCCACGCCCACTCCCCCAAGCGGGGAAGGAGCTTCGCGCTGCCGCTTGGCTGGGGACTGGGCACCGCCCTCCCGCGGCTCCTGAGCCGGCTGCCACCAGGGGGCGCGCCAGCGGTGTCCGGGAGCCTAGCGGCGCGTGTGCAGCGGCCAGTGCACCTGCTCTGGCCCTCGCCGCGGTCTCTGCCAGGACCCCGACGCCCAGCCTGACCCTGCCATTCAGCGGGGCTGCGGCTCCA GCCCGCCCGCCCG 4 4 1.00 0.0149855 0.949846 0 0.75 0.25 0

Though it may be difficult to tell from the above example, the desired output is a 17 column space-delimited file formatted like so:
  chrom startCoordinate endCoordinate repeatMotif numOfRepeatsInRef 1000refBasesBeforeStart 1000refBasesAfterEnd repeatSeqFromRef minFlankLeft minFlankRight repeatPurity defaultSlippage defaultStutter fractionAinMotif fractionCinMotif fractionGinMotif fractionTinMotif

Since I am working with non-human genomes, I cannot simply use the author's sample marker info files. I could not find explicit instructions for the creation of this markerInfoFile in popSTR's full journal article (https://academic.oup.com/bioinformatics/article/33/24/4041/2525679), so I emailed the authors. they instructed me to use Tandem Repeat Finder (in the Linux command prompt) with the following options:
./trf409.linux64 reference_genome.fasta 2 7 7 80 10 22 1000 -d -h -ngs > reference_genome_TRFmap

Here is the documentation for TRF:
https://blaxter-lab-documentation.readthedocs.io/en/latest/trf.html
I applied TRF to my reference genome, but the output was in the following format:
Chrom
    startCoordinate endCoordinate PeriodSize    numOfRepeatsInRef   SizeofConsensusPattern  PercentofMatchesBetweenAdjacentCopies   PercentofIndelsBetweenAdjacentCopies OverallAlignmentScore  fractionAinMotif fractionCinMotif fractionGinMotif fractionTinMotif Entropy repeatMotif repeatSeqFromRef    1000refBasesBefore  1000refBasesAfter

Here is an example of the first two lines of my current output file:
@ScyDAA6_1;HRSCAF=23
108 120 5 2.6 5 100 0 26 0 0 15 84 0.62 TTTGT TTTGTTTTGTTTT CAGTAAAGTCTTTCTTTCCTCTAACATAGAAAGTACTACTAGATTAGTGC TCTGTGTATGCTCTCTATTCTCAACCTCCAGATGCCCGTTCACACTGAGC
414 424 5 2.2 5 100 0 22 63 0 18 18 1.31 ATGAA ATGAAATGAAA TGAACTGGAATGTACATGATTGAATTTAAATTACTTCTTTAAAAAATTCC AGTATTGTGAATTGGTGCTAAATAAATAAACTGAATGAAAAATAACTCAC

It looks like my current TRF output is sorted with the chromosome names at the top, and the microsatellites found at each chromosome following underneath the name. I think I can include the appropriate chromosome names into each row with no problem, and I can easily rearrange the data into the appropriate columns. My main issue is with the information itself. As you can see, my current output includes data I don't need: period size, size of consensus pattern, percent of matches between adjacent copies, percent of indels between adjacent copies, overall alignment score, and entropy are not required for my desired output (unless some of these terms are synonymous with information included in the desired output and I am simply not aware). Furthermore, my output needs to include the minimum flanking on the left and right, repeat purity, default slippage, and default stutter. Which TRF options do I need to use for this information to be included in my output?


text:
Converting FASTA Sequences to SNP FASTA Format
I have two FASTA sequences. Each is from a different species' mitochondrial genome. The sequences are pasted below.
>speciesA_control_region_partial_sequence_mitochondrial
TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA
>speciesB_control_region_partial_sequence_mitochondrial
TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA

My goal is to extract all SNPs in these sequences and format them in SNP FASTA format, as shown below:
>rs11540652
GAACC[A/G]GAGGCCC

I specifically need this format (rather than, say, VCF) so I can search for RFLP sites using the following website:
http://bio.kuas.edu.tw/snp-rflping2/rflpUI.jsp
Does anyone know of an existing program which will accomplish this? Also, what additional materials do I need to detect these SNPs? I presume that I need assembled mitochondrial reference genomes for each species - anything else?

Do you have a VCF file? It should be possible to compute one using bcftools and it would give the position in the reference sequence so generating these “SNP” sequences could be done using this and the reference sequence. The full sequence for both samples will take longer to compute and is not necessary.

Your question is not very clear. Where does the flanking sequence in :
>rs11540652
GAACC[A/G]GAGGCCC

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the [] is from the reference... I don't know an existing existing program but here's a bespoke one.
speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

Gives:
>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.
I hope that helps.


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:
Samtools Index: Chromosome Blocks not Continuous
I am working with short-read whole-genome sequences from the NCBI's SRA. I have aligned and sorted all of my short-read sequences and am attempting to index each sequence into .bai format using samtools index, but am running into a couple of errors.
I unpacked the original .sra files in the following manner:
fastq-dump --defline-seq '@$sn[_$rn]/$ri' --split-files SRR6509138.sra

I then aligned each fastq paired-end read, removed duplicate reads, and converted each to bam format like so:
bwa mem xiphophorus_birchmanni_10x_12Sep2018_yDAA6.fasta SRR6509136_1.fastq SRR6509136_2.fastq | samblaster -e -r| samtools view -Sb - > blasted_SRR6509136.bam

Some of my files did not come as paired-end reads. For those files, I used the --ignoreUnmated flag on samblaster.
I then sorted each bam file like so:
samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

I am attempting to obtain a bai file for each bam using the following command:
samtools index sorted_SRR6649368.bam sorted_SRR6649368.bam.bai

This is the error I run into for the unpaired reads:
[E::hts_idx_push] Chromosome blocks not continuous
[E::sam_index] Read 'HWI-ST387:164:D0CJWACXX:4:1101:1129:13264' with ref_name='ScdB1pO_646;HRSCAF=880', ref_length=33092979, flags=16, pos=23260108 cannot be indexed
samtools index: failed to create index for "sorted_SRR791885.bam": No such file or directory

This is the error I run into for the paired reads:
[E::hts_idx_push] Unsorted positions on sequence #79: 11159020 followed by 11158717
samtools index: failed to create index for "sorted_SRR6649368.bam"

Can anyone help me figure out why this is happening?

Are there some instructions you are following that tell you to use these commands? Where has the idea that you need to use -n on the samtools sort command here come from?

samtools sort blasted_SRR6649368.bam -o sorted_SRR6649368.bam -n

These error messages indicate that the reads are not sorted by coordinate — in particular, that the reads mapped to ScdB1pO_646;HRSCAF=880 are not all together, and that the reads at positions 11159020 < 11158717 are not sorted by position on their chromosome.
This is because samtools sort -n has been used to sort the reads by name instead. Remove -n to sort by position, which is what is needed to prepare a BAM file for indexing with samtools index.

Thank you, it worked! As for why I used the -n flag, I am trying replicate the results of a journal article. The journal article did not provide explicit instructions for sorting/indexing the files, and this is my first time sorting/indexing files. I think "sort by reads" simply sounded like the correct thing to do, but it looks like I was incorrect


text:
Difference between paired-end, mate-pair and long read
I writing here because I have some questions for you.
I wondered what the essential differences were between paired-end, mate-pair and long read?
for me mate-pair and classic paired-end are both paired-end reads, with the difference that :

For classical paired-end: the insert size on classic paired-end is
smaller (about 500bp)
while the insert size of mate-pair is much longer (several Kb)
which allows to join the contiguous between them especially is it?

In revenge for the long-reads, I imagine that they are simply reads that are synthesized with a large read size but that do not allow like the maite-pair to make junctions between contigs?

They are all very different in separate regards, but they all refer to different wet-lab and sequencing protocols/technologies.
First, PE (paired end) reads are typically short (50-300) reads, most often Illumina HiSeq, MiSeq or NovaSeq protocols. Both pairs originate from a single fragment which is sequenced from either end:
------> PE1
-------------------------------- Fragment
                         <------ PE2
       ^^^^^^^^^^^^^^^^^^        Insert

In contrast mate pairs arise from a fragment that is circularized before sequencing. Now imagine the fragment is circularized:
                     -----> MP1
....--------------------------------.......
          MP2 <-----
                    ^Circularized (joined) here

Long reads typically refer to PacBio or Oxford nanopore data, with very long reads lengths (up to several hundred kb).

thank you very much for your time, as you can see ATpoint gave a link where I also asked a question in respond to a post, sorry for the duplication...


text:
Python: How to write duplicate sequences removed from fasta file to new file
I currently am using this code to remove duplicate sequences from the fasta file. However, I would also like to write a new file with only the removed duplicates as well as a count for how many times they appear. Please let me know if you have any suggestions, I am hoping to keep it written in one py file, but I am not opposed to running a second file to do so. I attended to write the appended strings to a new file but it looked slightly off...
from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   if str(record.seq) not in seen:
       seen.append(str(record.seq))
   records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")
```

Hi @User this has been nicely answered from different perspectives. Could you please send off one of the answers to this question as "accepted". The site has certainly helped here, so to reciprocate that help, this certainly helps.

Edit: user @RamRS has a more robust answer for actually using this. Please use that. Leaving answer up as OP appears to be interested in implementing as a learning activity.
Please let us know why this is not doing what you intend.
I haven't checked it, but something like this would work:
from Bio import SeqIO
import time
from collections import defaultdict

start = time.time()

seq_counts = defaultdict(0)
records = []

 for record in SeqIO.parse("b4r2.fasta", "fasta"):
   seq = str(record.seq)
   if seq in seq_counts:  # python is smart about dict keys
       records.append(record)

   seq_counts[seq] += 1

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

# print the counts
for seq in seq_counts:
    print(seq_counts[seq], seq)

print(f"Run time is {(end- start)/60}")

It will be a bit slow for large fasta files and I am sure that some people will be very upset that I'm using a defaultdict for this.
The code that you have would probably not work as you would want the records.append() call to be only run in the if statement testing for non-membership.
There is absolutely no reason to have to write more than one script for this.

There is absolutely no need to write even this script when tools like seqkit rmdup already exist.

Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.

Thank you! I've seen some other scripts that can do similar things but I wanted to make my own... :) Appreciate it a lot.

i am indeed upset that you use collections.defaultdict, would prefer collections.Counter ;) e.g. seq_counts = Counter(record.seq for record in SeqIO.parse(...

@RamRS I'll note that seqkit rmdup does not appear to fully satisfy the asker's question, which asked for a specific count of each sequence.

@Chris_Rands yeah totally, i just know the defaultdict syntax by heart and i'd have to go look up how to use a counter again...

@MaximilianPress -D and -d do that. If OP wants to save count against each sequence (as opposed to against each sequence identifier), the file would be unreliable as it would need to be a FASTA with an annotation added in the header. With the -D, OP should be able to use seqkit grep and seqkit duplicate to get their custom file.

@RamRS all right then, that is what they should do if they're using this in earnest. On the other hand, it does seem like they want to implement it themselves as an exercise, so I'll leave my answer up for now with an edit.

Agree with you 100%. This is great as a learning exercise (who among us did not write our own fasta deduplicator?), but these robust tools, as you rightly call them, are better tested so they catch edge cases that our code might not look out for.

Please do not reinvent the wheel - seqkit rmdup does what you want.
From the website:

rmdup
Usage
remove duplicated sequences by id/name/sequence

Usage:
  seqkit rmdup [flags]

Flags:
  -n, --by-name                by full name instead of just id
  -s, --by-seq                 by seq
  -D, --dup-num-file string    file to save number and list of duplicated seqs
  -d, --dup-seqs-file string   file to save duplicated seqs
  -h, --help                   help for rmdup
  -i, --ignore-case            ignore case

How does a -D file look like, you ask:
$ zcat hairpin.fa.gz \
    | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt

$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2

Hi @RamRS perfect thank you very much, that is really helpful


text:
dbnSNP frequency anomalies
Sometimes dbSNP reports very different allele frequencies for different large-scale genome projects e.g. between 1000 Genomes and GnomAD
rs11822440 1000Genomes A=0.4629    C=0.5371 GnomAD A=0.99997   T=0.00003
rs28531832 1000Genomes T=0.4379    C=0.5621 GnomAD T=0.99993   A=0.00007
rs78931618 1000Genomes A=0.7334    C=0.2666 GnomAD A=0.9999    T=0.0001
rs75716913 1000Genomes T=0.9006    C=0.0994 GnomAD T=0.99149   A=0.00851
Is this an issue with one or other of these projects data processing? Is something else going on here biologically or technically that might cause this?

do you know what populations were used to calculate the gnomad frequencies?

These are SNP frequencies - its not the same thing as whole genome GC content. What we are looking at here is selection pressure at given SNPs

@user438383 the pattern is consistently different between projects in global and sub-populations e.g. African, European

@Michael these are the same SNPs with different frequencies for different projects - but both on "normal human" so I would be surprised if there is that much selection via recruitment.

The other thing - have you double checked you are comparing the same genome build between the two (i.e. gnomad is gr38 and 1kg is gr37)?

allopatry .. it happens

I'll stick my neck out and simply state that allopatry alone would account for the pattern you are observing. Thus populations can exihibit variation in allelic profiles, particularly where Mendellian inheritance is operating (i.e. higher eukaryotes), singly due to geographic isolation.
This represents the Meyers view of population genetics and would account for variation within neutral theory, e.g. silent mutations. Obviously I'm guessing the populations were selected by geographic foci in generating the allelic frequencies.


text:
e-utils api: get chromosome position for each snp on a gene
So I have the following e-utils api-link:

https://eutils.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=gene&db=snp&id=5726&retmode=json

it returns all the SNP-ids for the gene with the id 5726 as an example. Additionally to the SNP-ids I also need the SNP position and the genomic distance. Also I would need the chromosome but I could query that from the gene as well rather than for each SNP. Do you know how I can add all these informations to the resulting json/xml?
The first one in the resulting list is the SNP with the id 1490381980. Further information about it can be found on dbSNP: https://www.ncbi.nlm.nih.gov/snp/?term=1490381980%5BUID%5D
I really want to avoid creating a request for each single SNP and would prefer to query everything at once. Due to the unavoidable performance-issues. The Tas2R38 doesn't have that many SNPs but there are other genes with way more.

Hi @cowboy_patrick could you plase state/clarify what you mean by "to query everything at once. Due to unavoidable performance issues.".

It would be great if I could get all the SNPs and their positions with a single api-call for a specified gene. Instead of creating a call for each SNP on the variation services or something similar

If your dataset is limited to human SNPs, and you know the list of SNPs, I would advise taking a look at the dedicated NCBI Variation Services.  This would be one request for each SNP, but if you fire off 1-2 per second, I think you'll be ok.
For example, for rs1490381980, you could try something like:
curl -X GET 'https://api.ncbi.nlm.nih.gov/variation/v0/refsnp/1490381980' \
   | jq '.primary_snapshot_data.placements_with_allele[] | "\(.seq_id) \(.alleles[0].allele.spdi.position)"'

It would be better to query for the SNP and parse the JSON in python, but jq helps provide an idea of what is possible.
Alternatively, you could download the JSON or VCF data from the NCBI dbSNP FTP site, and process out the information you need.