text:
What is the correct way to map Hi-C data with bwa mem?
Library Prep

I have a Hi-C library prepped using an enzyme that cuts at GATC, so it leaves GATCGATC as the junction sequences. This library was sequenced on a 2x150 PE Illumina run.

Data Pre-processing

The reads were adapter and quality trimming with trimmomatic paired-end mode, keeping reads that overlap.
Finally, the forward and reverse reads were trimmed such that all bases 3' of the junction sequence (GATCGATC). For example, the read:

5' AAAATT**GATCGATC**AATTATGGATAA 3' becomes: 5' AAAATT**GATC 3'

In this process, sequences that don't contain the junction sequences were left alone. Trimming junction sequences in this way makes it such that each read pair maps to only two places in the genome. No three-way or four-way chimeric mappings are possible since everything 3' of the junction sequence was removed.

My Question

What is the correct way to now map these trimmed Hi-C reads to a genome using bwa mem?

Progress so far/things I have tried

In the past I have just tried mapping using bwa mem with no special parameters using this command: bwa mem ref.fa hi-c_f.fq.gz hi-c_r.fq.gz. One concern is that the computations that bwa mem performs on insert sizes in this mode will cause some mappings to be faulty or that read pairs with chromosome-length information will be removed from the alignment because they fall on the tail of the insert size distribution (non-Gaussian, but log decay).
Another option to consider is mapping each file, hi-c_f.fq.gz and hi-c_r.fq.gz, independently then adding mate information with samtools fixmate. However, I do not know how to do this.
A last option that I have considered is to use bwa mem in paired-end mode, but to use a special flag. In the past I believe that -5 was the correct flag, but in this post Heng Li discusses that the -5 flag behavior has changed, as well as a new -q flag. In this post Heng Li discusses that he thinks bwa mem is still a good option for mapping Hi-C data but does not discuss the proper flags.

Limitations

The answer should provide the correct way to map Hi-C using bwa mem, but other software suggestions will also be considered.

This ends up depending on the software you want to use downstream of mapping. For example, HiCExplorer still needs reads to be aligned as single-end. I imagine it's not alone in that regard.

See here, in particular slide #10 from the tutorial:
bwa mem -5SPM ref.fa read1.fq read2.fq > out.sam

Here, -SP disables pairing. The Aiden lab and this paper also use a similar command line. Beware that there are many Hi-C pipelines, but not many are using bwa-mem.

Thanks for your post. It looks like the -SP flags are sufficient to map the trimmed data as I discussed above, while -5SPM is appropriate for Hi-C data that do not have the linker sequence trimmed.

There are some interesting suggestions from the HiCExplorer developers:

Mates have to be mapped individually to avoid mapper specific
  heuristics designed for standard paired-end libraries.

# map the reads, each mate individually using
# for example bwa
#
# bwa mem mapping options:
#       -A INT        score for a sequence match, which scales options -TdBOELU unless overridden [1]
#       -B INT        penalty for a mismatch [4]
#       -O INT[,INT]  gap open penalties for deletions and insertions [6,6]
#       -E INT[,INT]  gap extension penalty; a gap of size k cost '{-O} + {-E}*k' [1,1] # this is set very high to avoid gaps
#                                  at restriction sites. Setting the gap extension penalty high, produces better results as
#                                  the sequences left and right of a restriction site are mapped independently.
#       -L INT[,INT]  penalty for 5'- and 3'-end clipping [5,5] # this is set to no penalty.

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R1.fastq.gz 2>>mate_R1.log | samtools view -Shb - > mate_R1.bam

$ bwa mem -A1 -B4  -E50 -L0  index_path \
    -U mate_R2.fastq.gz 2>>mate_R2.log | samtools view -Shb - > mate_R2.bam

# build matrix from independently mated read pairs
# the restriction sequence GATC is recognized by the DpnII restriction enzyme

$ hicBuildMatrix --samFiles mate_R1.bam mate_R2.bam \
                 --binSize 10000 \
                 --restrictionSequence GATC \
                 --threads 4
                 --inputBufferSize 100000
                 --outBam hic.bam \
                 -o hic_matrix.h5
                 --QCfolder ./hicQC

Option -SP exactly aims to disable "mapper specific heuristics designed for standard paired-end libraries".


text:
How to align genomic sequence with corresponding amino acid sequence
Does anyone know of a program that can align a genomic sequence with introns with the corresponding amino acid sequence?
I have both the genomic sequence and the correct amino acid sequence but no information on the genemodel in e.g. gff or genbank format.
This is not what one would call a 'normal' alignment, but I would like it to look like this:
Hypothetical example:
nucleotide sequence with intron in italic:
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
corresponding amino acid sequence:
MHDTTDVRTYVRRRRL
The intron does not encode for amino acid which should thus result in a gap in the 'alignment'.
Desired 'alignment':
ATGCACGATACGACTGACGTACGTACGTACGTACGTACGTACGTACGGAGACGTAGACTC
  M  H  D  T  T  D  V  R              T  Y  V  R  R  R  R L

Have a look at exonerate's protein2genome!
From the documentation:

Aligning a protein to genomic sequence:
Similarly, it is possible to align a protein sequence to the genome, (similar to GeneWise, but with
  heuristics).
exonerate --model protein2genome query.fasta target.fasta
This model will allow introns in the alignment, but also allow frameshifts, and
  exon phase changes when a codon is split by an intron.

The classic tools for this are genewise and the already mentioned exonerate. Both are capable fo aligning protein sequences to genomic DNA (or cDNA) while modelling splice sites to give you a gene model.
For example, this is the output1 of genewise on your sequences:
EMBOSS_001         1 MHDTTDVRTYVR----RRRL
                     MHDTTDVRTYVR    RRRL
                     MHDTTDVRTYVRTYVRRRRL
EMBOSS_001         1 acgaaggcatgcatgcacac
                     taaccatgcatgcatggggt
                     gctgtcatgcatgcagatac

The output also includes a GFF of your sequence:
EMBOSS_001  GeneWise    match   1   60  22.88   +   .   EMBOSS_001-genewise-prediction-1
EMBOSS_001  GeneWise    cds 1   60  0.00    +   0   EMBOSS_001-genewise-prediction-1

Back when I was using these tools daily to align proteins to genomes, I found that exonerate was much faster but sometimes slightly less accurate. This was a few years ago now though, so that might have changed. I would always suggest trying exonerate first and only going for genewise if it fails.

1This link will only be valid for a short while, but I just used the sequences from the OP as input for the GeneWise webserver.

Thanks, The GeneWise seems to work as well, but the Resulting GFF is quite useless. It would be great if it would give an annotation for the intron, but that is unfortunately not the case.

@user1817 what do you mean? Of course it annotates the intron. It won't in your example since that's too small to be a real intron, but try it on a real sequence and the GFF will include the introns. I just ran p54 against its gene, for example, and the introns are there.

source: Stack Exchange

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

sort: votes


text:
Inspection of gene expression in scRNA-seq data
I am running the data preprocessing pipeline for scRNA-seq data presented here.

3.8.6.1 Gene expression
In addition to removing cells with poor quality, it is usually a good idea to exclude genes where we suspect that technical artefacts may have skewed the results. Moreover, inspection of the gene expression profiles may provide insights about how the experimental procedures could be improved.
It is often instructive to consider the number of reads consumed by the top 50 expressed genes.

By analyzing my data I have the following most expressed 50 genes (evaluating the total number of transcripts across cells):
Original

Log-transformed

Since I'll be working on log-transformed data, I assumed that the distributions of the top 50 genes were not excessively skewed. However, as you can see in the first plot, there are some samples (cells) with very high number of transcripts for the top gene. Is it something I should worry about? What other insights can I get from these two figures?
Edit
Original analysis (figure attached) was:

The distributions are relatively flat indicating (but not guaranteeing!) good coverage of the full transcriptome of these cells. However, there are several spike-ins in the top 15 genes which suggests a greater dilution of the spike-ins may be preferrable if the experiment is to be repeated.

In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the point the course instructors are trying to get across with that message box.
Eyeballing the results, you should also exclude any ERCC genes as well, as they are synthetic constructs and not from the target samples (although record their counts, as it may be informative for QC/normalisation purposes); they'll have odd distributions because of this.
ERCC sequences are typically added based on the concentration of DNA in a sample, rather than the number of cells. However with single-cell data (i.e. with known cell counts), it'd be possible to do a proper ERCC-based normalisation, assuming an equal amount were added for each cell. I suspect that knowing a bit more about the ERCC addition protocol would help in understanding why their distribution is so different.

Thanks! I noticed that the original comment for the analysis was not included in the question: I updated the question, I didn't get what they meant with flat distribution.

I forgot to ask, why remove ribosomal genes? I tried to find an explanation of this contamination, but with no results.

As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite difficult to completely remove from a sample (even with polyA selection), with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to. This high, variable expression throws out normalisation calculations, so it's better to remove them first.

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


text:
How to subset a BAM by a list of QNAMEs?
I have a text file 'qnames.txt' with QNAMEs in the following format:
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

I would like to subset my BAM file.bam via all of these QNAMEs into a new SAM.
Naturally, I can do this individually, e.g.
samtools view file.bam | grep 'EXAMPLE:QNAME1' > subset.bam

But I'm unsure how to do this for a list of QNAMES:

How do I write a for loop that will do all of these queries, outputting the correct SAM needed?
I could write a for loop that creates n SAM files, and then cats them…

Is there a way to specifically grep by QNAME? The above may grep reads that may not be associated with the correct QNAME?

How do I keep the BAM header?

Are these QNAMEs completely unrelated? If, by chance, you simply want to select all reads that belong to the same read group (which is determined by the QNAME), then there’s a much more efficient and simpler solution.

@KonradRudolph mind adding that solution over here?

@Beeba Given that OP hasn’t confirmed whether splitting by read groups is indeed desired, I’m reluctant to post an unrelated answer here. That said, it’s simply what the command samtools split does.

@KonradRudolph makes sense, no worries. I will look into samtools split. Thanks

(1)
It won't be super fast but you can provide grep with a file of QNAMES.
samtools view file.bam | grep -f 'qnames.txt > subset.sam

where qnames.txt has
EXAMPLE:QNAME1
EXAMPLE:QNAME2
EXAMPLE:QNAME3
EXAMPLE:QNAME4
EXAMPLE:QNAME5

(2)
This would be a little more complicated but can you give an example where the grep might be have the correct QNAME?
(3)
To keep the BAM header I would use samtools -H file.bam > header.txt to get the header and then cat the header file with the grep'ed file

This is quite slow if the bam file or the list of reads to keep are of appreciable size.

Since the qnames are all fixed strings, using fgrep instead of grep will provide a significant speedup — several orders of magnitude. However, both approaches may lead to spurious results: If you have a qname foo in our list, and the BAM file also contains qnames fooX and fooY, the latter will also be returned in your results, erroneously.

As it was already said, this is very slow. You can easily do this much faster with samtools or picard (see mine or Pierre's answers).

use picard FilterSamReads http://broadinstitute.github.io/picard/command-line-overview.html#FilterSamReads

READ_LIST_FILE (File) Read List File containing reads that will be included or excluded from the OUTPUT SAM or BAM file. Default value: null.

UPDATE 2021/06/28: since version 1.12, samtools now accepts option -N, which takes a file containing read names of interest.
Using samtools 1.12 or greater:
samtools view -N qnames_list.txt -o filtered_output.bam input.bam

OLD ANSWER:
When it comes to filter by a list, this is my favourite (much faster than grep):
# given a bam file "aln.bam" and a list of read names "reads.txt":
samtools view -h aln.bam \
| awk 'FNR==NR {reads[$1];next} /^@/||($1 in reads)' reads.txt - \
| samtools view -b - > filtered_aln.bam

Basically, the awk command first parses the reads.txt file, creates an array of QNAMEs, then parses the alignments and print if the QNAME is present in the array (or if the line belongs to the header).

Below is a small bit of python code showing one naive but manageable way of doing this (N.B., I expect grep to be faster, though getting it to output the header will be annoying):
import pysam
qnames = set(...)  # read names go here
bam = pysam.AlignmentFile("some input file.bam")
obam = pysam.AlignmentFile("some output file.bam", "w", template=bam)

for b in bam.fetch(until_eof=True):
    if b.query_name in qnames:
        obam.write(b)
bam.close()
obam.close()

Now that'll work, but if you have a lot of names to look through then it'll end up being quite slow (as will grep -f ...). If you need to select a large list of reads by name, then a more efficient strategy is:

Name sort the BAM file (samtools sort -n or picard), which can be done with multiple threads.
Sort the list of read names to match (this will be more complicated than one would naively think, since samtools and picard will name sort differently, so be sure to put some thought into this).
Perform the same sort of iteration as above, but only compare to the top read name in your list from step 2 (removing that element from the stack after a match).

Keeping the BAM header with the grep approach is simply a matter of echoing it first: (samtools -H input.bam; samtools input.bam | grep …) | samtools -b - -o output.bam

Performing a grep on $n$ alignments to $m$ qnames would give you $O(mn)$ operations. If you decide to sort your BAM and qnames first, you can reduce it to $O(\min\{m, n\})$ by popping reads and qnames off their stacks. You also don't read the entirety of the BAM at once this way, limiting memory consumption. I'm pretty sure Picard does something like this, as well.
Here's an example in Python that I've written for this very purpose:
import pysam

def filter_qname(bamfile, idfile, outfile):
    '''
    Filter reads in a SAM/BAM file by their query names

    Parameters
    ----------
    bamfile : str
        Sorted BAM file path
    idfile : str
        Text file path containing qnames to keep
    outfile : str
        Output file to write to
    '''
    # read in name-sorted BAM file
    bam = pysam.AlignmentFile(bamfile, "rb")
    if outfile.endswith('.sam'):
        output = pysam.AlignmentFile(outfile, "w", template=bam)
    elif outfile.endswith('.bam'):
        output = pysam.AlignmentFile(outfile, "wb", template=bam)
    else:
        raise ValueError("Unknown output file format for `outfile`: {}".format(outfile))
    #  read in IDs to be removed (list(set(...)) to only keep unique IDs)
    ids = list(set([l.rstrip() for l in open(idfile, 'r').readlines()]))
    # sort IDs to match BAM for efficient processing (destructive function)
    ids.sort(key=str.lower, reverse=True)
    n_ids = len(ids)
    # variable for ensuring BAM is sorted by query name
    last_q = None
    reads = bam.fetch(until_eof=True)
    read = next(reads)
    while True:
        # if read name is greater than current top of stack
        if read.query_name > ids[-1]:
            # if this is the last ID
            if n_ids == 1:
                # write remaining reads to new file and break out of while loop
                output.write(read)
                for r in bam:
                    output.write(read)
                break
            # otherwise pop that ID, try again
            else:
                ids.pop()
                n_ids -= 1
        # skip reads that match the top of the stack
        elif read.query_name == ids[-1]:
            read = next(reads)
        # if read name is less that top of stack, write it and move on
        else:
            output.write(read)
            read = next(reads)
    output.close()

I've implemented it in my own miscellaneous bioinformatics tools package here if you want to see it in action

Most of the solutions posted thus far are either slow, or can yield spurious results in certain cases. (I haven't tested Picard, which I assume will work as intended. But personally I tend to recoil when I have to set up a JVM!)
Let's approach the problem from a different angle:
Good hash functions have ~O(1) complexity for a lookup. If we can create a hash table of our QNAMEs, then searching a BAM file should be O(n) in the number of records. The Python developers are very proud of the hash function in cpython, so let's use that:
#!/usr/bin/env python3

import sys

with open(sys.argv[1], 'r') as indexfile:
    ids = set(l.rstrip('\r\n') for l in indexfile)

for line in sys.stdin:
    qname, _ = line.split('\t', 1)
    if qname in ids:
        sys.stdout.write(line)

With a test set consisting of a ~3GB BAM file and a qnames.txt with ~65K entries, running samtools view input.bam | ./idfilter.py qnames.txt takes roughly 6.5 seconds on my test server. By comparison, piping the SAM to fgrep -f qnames.txt (using GNU grep) takes about 8.5 seconds (but may yield spurious results).
(Why might grep yield spurious results? If you have a QNAME foo in your qnames.txt, but your BAM file also contains fooX and fooY, they will all be matched by fgrep. The solution is to anchor each of your search patterns between the beginning of the line and a tab, e.g. ^foo\t, but then you have to use standard grep which will have to construct an NFA for each pattern, and not fgrep which implements Aho-Corasick, and you'll be orders of magnitude slower.)
I've rewritten my little Python script in Rust using std::collections::HashMap for the hash function, and using the rust_htslib crate to read and write straight from BAM, and this is more than twice as fast as the Python, even when reading and writing BAM. However, my Rust isn't really fit for public consumption…

This is a great answer! It's quite comprehensive, and I agree with your distaste of setting up a JVM :)

You can use SAMsift:
samsift -i file.bam -0 'q=open("qnames.txt").read().splitlines()' -f 'QNAME in q'

-i file.bam specifies the input BAM file. -0 'q=open("qnames.txt").read().splitlines()' loads your list of QNAMEs. -f 'QNAME in q' specifies the filter for alignments.
With SAMsift one can use any Python expression for filtering. The downside is that it's slower than the other tools.

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


text:
Can I run STAR without an annotation file?
I wish to use Rascaf to scaffold a fragmented draft genome.
For this, I need to provide a BAM file of aligned RNA-seq reads and the draft genome.
So, I indexed the draft genome with STAR like this:
STAR --runMode genomeGenerate --genomeDir output/index --genomeFastaFiles draft.fa

Then, I tried to aligned the reads like this:
STAR --genomeDir output/index --readFilesIn reads.fastq.gz --readFilesCommand gunzip -c --outFileNamePrefix output/alignment --quantMode GeneCounts --outSAMunmapped None --outSAMtype BAM SortedByCoordinate --outSAMattrRGline ID:RG1 CN:yy \"DS: z z z\" SM:sample

The latter commands outputs the following error:
Transcriptome.cpp:51:Transcriptome: exiting because of *INPUT FILE* error: could not open input file exonGeTrInfo.tab
Solution: check that the file exists and you have read permission for this file
          SOLUTION: utilize --sjdbGTFfile /path/to/annotantions.gtf option at the genome generation step or mapping step

It suggests I should use a GTF file. But this is a draft genome of an individual for which no GTF is available. I could try to adapt a GTF from a close relative but for the scaffolding, it seems superfluous. Is there a way I can map with STAR without the GTF?
These commands were run on a GNU/Linux machine.

I don't think you can use the --quantMode GeneCounts option with no annotations. I think the error is trying to look for an exon file generated from the annotations to do the quantitation on. Remove that and I think it should work, as the manual specifically states that annotations are optional but highly recommended.

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


text:
How do I rename fasta headers?
I am trying to use the seqkit replace command to replace chromosome names in the format chr_I, chr_II, ... to I, II, .... I am using the following command:
seqkit replace -p "(.)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

My keyValues.txt file contains the following:
chr_I   I
chr_II  II
...

The two columns are separated by a tabulation.
I get the following output:
[INFO] read key-value file: keyValues.txt
[INFO] 6 pairs of key-value loaded
[ERRO] pattern "(.)" matches multiple targets in "chr_I", this will cause chaos

I don't know what I am doing wrong. I chose the (.) pattern to match the whole header but it seems to be wrong. Any help would be appreciated.
Update
All the headers are shown below:
grep '>' mySequence.fasta
>chr_I
>chr_II
>chr_III
>chr_IV
>chr_V
>chr_X

Show some of the headers, the pattern argument is probably not suitable. I am guessing that you'll get headers with duplicate names now (which will cause chaos).

I have written all the header names in my question.

There are 4 matches to your pattern "chr_I"

Oh, so the key is a regular expression? I tried putting a dollar at the end of the key names, but I get the same error.

Why not use sed instead?
sed -e 's/chr_I/I/' -e 's/chr_V/V/' -e 's/chr_X/X/' mySequence.fasta > mySeq.fasta

Or even simpler:
sed 's/chr_//' mySequence.fasta > mySeq.fasta

Worked. I just don't know sed. I need to learn it.

Great! It works a bit like grep (which you know how to use), but then with replacement function.

I would suggest using sed 's/^chr_//' instead. Anchoring the pattern to the beginning of the line should speed things up for large files. Otherwise, sed will have to scan the entire length of each line in case there's a match.

@terdon, good point. But does it work with fasta? Remember the > at the beginning of the header.

@b.nota Duh! I should learn to read :) Indeed, sed 's/^>chr_/>/ would be better.

Minor correction for the last: sed 's/^>chr_/>/' ( in order to keep the >)

I like using awk when there's any kind of lookup involved: something like
awk 'FNR==NR {
    hash[">"$1]=$2;
    next;
}

if ($1 ~ /^>/) {
    print ">"hash[$1];
} else {
    print $1;
}' keyValues.txt mySequence.fasta

Basically the FNR==NR check tells awk to work only on the first file (useful primer here), keyValues.txt and create an association for each key (eg. chr_I) with its value (eg. I). The rest of the code after the next works only on mySequence.fasta, printing out the lookup value only if the line is a fasta header, as checked by the $1 ~ /^>/ condition.

Your issue is that the pattern . matches a single character. To match the whole header, use .*:
seqkit replace -p "(.*)" --replacement "{kv}" --kv-file keyValues.txt mySequence.fasta

Alternatively as pointed out by others you may not really need keyValues.txt here as you're just removing the chr_ prefix:
seqkit replace -p "chr_(.*)" --replacement '$1' mySequence.fasta

(Don't forget the single quotes for $1 else the shell will expand it)

Here is an alternative to sed: cat test | tr -d "chr_"
The pipe is to output the file, while the tr command is to translate and change the file

The use of the cat command is completely redundant here (and imposes a runtime overhead). Just do tr -d 'chr_' < test.


text:
Duplicate gene symbol handling in GEO gene expression data
I have downloaded a gene expression data from GEO database (GSE3268) in which in some of its rows there are duplicate gene symbols.
For example TP53 exists in two rows with different expression values and different GenBank Accession IDs.
They have the following Nucleotide Title:

Human p53 cellular tumor antigen mRNA, complete cds
Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA

I want to map these data with edge list of gene symbols to make a network. Which expression value should I consider for such duplicate gene symbols for mapping?

Since this is a microarray (U133A), what you're seeing is that each probe is associated with one or more transcripts. The general strategy when dealing with microarrays is to use RMA (or fRMA) to summarize to probesets and then collapseRows to summarize to gene-level information. Since this is all typically done in R, you can find some further discussion on the bioconductor site. I would suggest going through this process and then only using the result of that in building your network (whether you're happy with sticking to transcript-level quantifications or gene-level quantifications for building your network is up to you).

I have a problem with using fRMA since the there is no supplementary .cel file in GEO database for this dataset. How can I use R on .SOFT files to summarize gene information?  Another point is that in one of the subsets of my data set, (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM73388), it seems that RMA have been used and data is collected in .SOFT file, but again there is duplicates in it.

If RMA has already been used then you needn't run it a second time.

Very useful answer! In EuropePMC I found 63 publications mentioning the use of collapseRows(). Many  don't report the setting they apply. Some high impact publications reported as follows: "MaxMean" (default) = Gandal (Science2018), Raznahan (PNAS2018), Langfelder (NatNeurosci2016); "maxRowVariance" = Langfelder (NatNeurosci2016)


text:
Check if a set of genes fall more on a specific cluster than other clusters
I have a network which I clustered into 11 smaller clusters. I have a set of genes and I would like to check if statistically significant number of these genes fall on any specific cluster. One method which I thought of was to find which percentage of a total cluster were made up of my genes of interest. For example, say 60% of the genes in a cluster were the genes in my list. However, I was hoping that there was a better way to do this. This study on honey bees has done something of the kind.
I would appreciate if there was a way to do something similar. I do not know how to replicate what they have done in the study. Any help is appreciated.

Do you have problem with the hypergeometric test or with the later methods they do?

@Llopis I have never done hypergenometric tests before and I am not sure how to do it.

A hypergeometric test is related to the fisher test, in a table with one  dimension being some factors and the other one some other factors it is tested if they are independent. In R you con do this:
data <- matrix(c(0, 1, 2, 5), ncol = 2, nro = 2)
fisher.test(data)

Or you can use the specific function phyper, to see the relation between them you can read this question and answer


text:
Visualize similar interactions in two homologous networks
Background
I have two networks with homologous genes from two different species.
To identify homologs I performed a functional comparison. One of the species is not well studied while the other one is the house fly (D. melanogaster). I used ClueGO to predict the functions of the genes in the D. melanogaster network.
I found a function of interest influenced by a set of genes and predicted that the same genes in the other network (thus far uncharacterized) should also perform a similar function.
Questions
Q1: How do I check if there are overlapping patterns between both the networks on Cytoscape?
Q2: There are more than 100 homologs in this list and I would like to minimize the list by only considering those homologous genes that are arranged the same as well.

Do you hope to perform any more analyses on these networks aside from visualization?

@conchoecia I have edited the question. Please check.

Thanks for your edits. The part that states, "minimize the list by only considering those homologous genes that are arranged the same as well", is still unclear. It is clear that you want to narrow the list of homologs but it is not clear what "arranged the same" means. Does this mean that they have the same synteny in the genome, the same edges and edge weights in a network, or something else? We can also provide you with a better answer if you post text examples of your data.

Possible duplicate of Compare two networks of the same genes between two species

@gabrielet I did ask that previous question as well. That was about inferences that can be drawn from topological comparisons.

@conchoecia The edges have no weights but I would say consider those homologues with similar edges would more than suffice.

@TheLastWord, I see the difference but I believe that the goal is the same in both questions since what you want to achieve is the comparison of networks to find similarities.

An easy way to visualize this is to make a tab-delimited table that contains edge information for both species and both networks, color the edges depending on the species, and make the thickness and transparency dependent on the interaction strength.
For example, say you have a list of orthologs for the two species and just assign them a number or a name. Each edge has a "source" and a "target". The edge weight is the interaction strength.
(Note do not copy-paste the below table as the mixed tabs/spaces may make the graph topology incorrect). In my toy dataset the data look like this:
source  target  variable    value
gene1   gene2   speciesA    3
gene1   gene3   speciesA    42
gene2   gene3   speciesA    3
gene1   gene2   speciesB    9
gene2   gene3   speciesB    5
gene2   gene4   speciesB    41
gene3   gene4   speciesB    32

Import your network to cytoscape as a network from a file:

Then go to the style tab:

... and make sure edge is selected.

Then you can adjust the transparency, edge color, and edge stroke settings, like so using edge stroke color:
:
... and transparency panels as examples.

Don't forget to save your style when it is complete!
When you have adjusted these settings you will have something that looks like this:

I think what you're looking for, as I suggested in this old question (quite similar to this one, I believe), is the GASOLINE app for Cytoscape that allow you to perform an analysis that will tell you which parts of the two networks are conserved and which genes belong to the overlapping subnetworks.


text:
cell type specific genes - heatmap using rank based approach
I was wondering if some one can help me with pseudo r code for my approach here based on this paper to plot a heatmap
Figure 3

Figure description
The top 40 enriched genes percelltypeare shown in a heat map.Onlyhighlyexpressed genes with FPKM ⬎20 areincludedin this analysis. Fold enrichment iscalculatedasFPKM of one cell type divided by the average FPKM of all other cell types. The majority of these genes showed specific expression by only one cell type, with the exception that some are expressed during more than one maturation stage in the oligodendrocyte lineage.
Similar description in their methods:
Differential expression was calculated as the FPKM of a given cell type
divided by the average FPKM of all other cell types. Genes were ranked by
their fold enrichment in each cell type, and top enriched genes for each
cell type were identified
In my case I have 9 cell types(columns) and ~24000 genes(rows) from mouse and the values are TPMs.
I know to order celltype only based on one column like this
expr.order = expr[order(expr[,1],decreasing =T),]

Which clearly not the right way and would need some help with code in r

It’s unclear from your text what you actually want to do: there’s no question. For ordering the expression values based on one column you already found the solution. If that’s not what you mean you need to be more specific.

I do not know how to achieve this Genes were ranked by their fold enrichment in each cell type  which I guess means rank individual column independently not sure how to do this

You can't rank individual columns independently (of the same data), but you can order them by several factors.

Well in that case you first need to compute fold changes using one of the existing differential expression packages, e.g. DESeq2 or edgeR.

BTW, what else have you tried with R? Did you compute the fold change and found the top 40 as described in the image?

The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:
> expr <- matrix(runif(9*24000), nrow=24000, ncol=9)
> top_genes_each_celltype <- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
> dim(top_genes_each_celltype)
[1] 40 9
> top_genes_all_celltypes <- unique(as.vector(top_genes_each_celltype))
> length(top_genes_all_celltypes)
[1] 358

Then you can build a biclustered heatmap with expr[top_genes_all_celltypes,], for example, using heatmap2.


text:
What are negative FSC-A values in flow cytometry data?
I did flow cytometry analysis of E. coli cells 7 hours post inoculation (1% overnight culture in 3ml).
I am analysing the FCS3.0 files using the flowCore package in R. I observe that FSC-A has negative values in some events, although the corresponding FSC-H values were all positive. These points comprised just a minor fraction of the data (~250/10000 in this case). The FSC voltage was set at 788V and the threshold flow rate was 5000-10000 events/sec.
What do these negative values mean?

FSC-A vs FSC-H: plotted only negative FSC-A values

SSC-H vs FSC-A: red points indicate negative FSC-A values

Why did you delete the SSC-H plot? Also, yould you post the corresponding SSC-A/-W FSC-W plots?

@bukwyrm well, what I claimed in the image that I posted wasn't really true. I don't have the -W data. I tried to see if there is any correlation but the only thing that I could find is that these negative FSC-A points have low SSC-H  but removing these points from the total data does not change the distribution of the latter.

This is probably a matter of calibration. If this is not the raw data (actual voltage curves), look at that.
As an area below 0 is not possible, it has to be a calibration issue, either in the machine, or in the software.
the -H plots are about the maximum height of the voltage spike in the side or forward light detector, while the -W plots are about the width (using some definition of flank-detection, so variation prone).
The area is also computed in various ways: using something like (H*W)/c (with c something like 2) or with an actual computation of the area of the spike, but again this has all kinds of predefinitions in it (what's the beginning/end of a  spike, what's the background-level, ...) Some of these predefinitions you can alter in software, some are hardcoded or even part of the functionality of the hardware, so it might be easy to impossible to get rid of the negative values.
The FSC-H values corresponding to the FSC-A values are in your plot, and positive, so the most obvious solution (a bad algorithm that calculated Area as Height times Width times Magic and did not account for negative Heights (which depending on spike definition, you might get)) is out.

So basically there is some unknown bad algorithm that is responsible for this. Should I consider FSC-H as the reliable estimate and retain these points for further analysis or discard them? As I mentioned in the other comment, removing these points does not change the SSC or FITC distribution.


text:
Understanding DESeq2 design, contrast and results
I have a set of high-troughput experiments with 2 genotypes ("WT" and "prg1") and 3 treatments ("RT", "HS30" and "HS30RT120"), and there are 2 replicates for each of the genotype x treatment combinations.
The read counts for the genes are summarized in a file that I load as follows in R:
> counts_data <- read.table("path/to/my/file", header=TRUE, row.names="gene")
> colnames(counts_data)
 [1] "WT_RT_1"          "WT_HS30_1"        "WT_HS30RT120_1"   "prg1_RT_1"
 [5] "prg1_HS30_1"      "prg1_HS30RT120_1" "WT_RT_2"          "WT_HS30_2"
 [9] "WT_HS30RT120_2"   "prg1_RT_2"        "prg1_HS30_2"      "prg1_HS30RT120_2"

I describe the experiments as follows:
> col_data <- DataFrame(
    geno = c(rep("WT", times=3), rep("prg1", times=3), rep("WT", times=3), rep("prg1", times=3)),
    treat = rep(c("RT", "HS30", "HS30RT120"), times=4),
    rep = c(rep("1", times=6), rep("2", times=6)),
    row.names = colnames(counts_data))
> col_data
DataFrame with 12 rows and 3 columns
                        geno       treat         rep
                 <character> <character> <character>
WT_RT_1                   WT          RT           1
WT_HS30_1                 WT        HS30           1
WT_HS30RT120_1            WT   HS30RT120           1
prg1_RT_1               prg1          RT           1
prg1_HS30_1             prg1        HS30           1
...                      ...         ...         ...
WT_HS30_2                 WT        HS30           2
WT_HS30RT120_2            WT   HS30RT120           2
prg1_RT_2               prg1          RT           2
prg1_HS30_2             prg1        HS30           2
prg1_HS30RT120_2        prg1   HS30RT120           2

I want to build a DESeq2 object that I could use to either:

find differentially expressed genes when the treatment varies for a given fixed genotype

or:

find differentially expressed genes when the genotype varies for a given fixed treatment

In the bioconductor help forum I think I've found a somewhat similar situation, and I read the following:

Try a design of ~ genotype + genotype:condition
Then you will have a condition effect for each level of genotype, including the reference level.
You can constrast pairs of them using the list style of the 'contrast' argument.

However, this doesn't explain how to apply this "list style" to the "contrast" argument. And the above situation seems to be asymmetrical. By that I mean that genotype and condition do not seem to have an interchangeable role.
So I tried the following more symmetric formula:
> dds <- DESeqDataSetFromMatrix(
                                countData = counts_data,
                                colData = col_data,
                                design = ~ geno + treat + geno:treat)
> dds <- DESeq(dds)

Now, can I for instance get the differential expression results when comparing treatment "HS30" against "RT" as a reference, in genotype "prg1"?
And how?
If I understand correctly, the above-mentioned "list style" uses names given by the resultsNames function. In my case, I have the following:
> resultsNames(dds)
[1] "Intercept"               "geno_WT_vs_prg1"
[3] "treat_HS30RT120_vs_HS30" "treat_RT_vs_HS30"
[5] "genoWT.treatHS30RT120"   "genoWT.treatRT"

I guess I would need a contrast between "genoprg1.treatRT" and a "genoprg1.treatHS30", but these are not in the above results names.
I'm lost.

It seems that the "combining factors" trick described in part 3.3 of DESeq2 current "vignette" (as of may 2017) under the title "Interaction" is a way to access to the desired contrasts.
It seems possible do do it directly when building the colData and when calling DESeqDataSetFromMatrix:
Let's add a combined "geno" and "treat" factors to the future colData parameter:
> col_data$geno_treat <- as.factor(paste(col_data$geno, col_data$treat, sep="_"))
> col_data
DataFrame with 12 rows and 4 columns
                        geno       treat         rep     geno_treat
                 <character> <character> <character>       <factor>
WT_RT_1                   WT          RT           1          WT_RT
WT_HS30_1                 WT        HS30           1        WT_HS30
WT_HS30RT120_1            WT   HS30RT120           1   WT_HS30RT120
prg1_RT_1               prg1          RT           1        prg1_RT
prg1_HS30_1             prg1        HS30           1      prg1_HS30
...                      ...         ...         ...            ...
WT_HS30_2                 WT        HS30           2        WT_HS30
WT_HS30RT120_2            WT   HS30RT120           2   WT_HS30RT120
prg1_RT_2               prg1          RT           2        prg1_RT
prg1_HS30_2             prg1        HS30           2      prg1_HS30
prg1_HS30RT120_2        prg1   HS30RT120           2 prg1_HS30RT120

We can now use a design where differential expression will be explained by these combined factors:
> dds <- DESeqDataSetFromMatrix(
    countData = counts_data,
    colData = col_data,
    design = ~ geno_treat)

We run the analysis:
> dds <- DESeq(dds)

Then we can query results for a particular contrast between such factor combinations. For instance, to have the results for the effect of treatment "HS30" against the reference state "RT" in genotype "prg1":
res <- results(dds, contrast=c("geno_treat", "prg1_HS30", "prg1_RT"))

The simplest manner is to not use a wald test, but rather an LRT with a reduced model lacking the factor of interest:
dds = DESeq(dds, test="LRT" reduced=~geno+geno:Treatment)

The above would give you results for Treatment regardless of level while still accounting for a possible interaction (i.e., a "main effect of treatment, regardless of the type of treatment").
As an aside, this is probably a case where the edgeR-preferred way of creating groups of genotype-treatment combinations and then using a model of ~0+group might make your life a bit easier. You'll get the same results (more or less) regardless, but it'll probably be easier for you to think in those terms rather than remembering that the base level will be treatment HS30 and geno prg1.

I tried the following:
dds <- DESeqDataSetFromMatrix(countData = counts_data, colData = col_data, design = ~ geno + treat + geno:treat) then dds = DESeq(dds, test="LRT", reduced=~geno+geno:treat) but this fails with Error in nbinomLRT(object, full = full, reduced = reduced, betaPrior = betaPrior,  : less than one degree of freedom, perhaps full and reduced models are not in the correct order. Besides, what I want is not "regardless of the type of treatment". I want results for a given fixed treatment or for a given fixed genotype. Maybe my question was not clear.

I'd have to go over the model matrix to see why it's producing the error. Anyway, I would suggest going the route the edgeR authors prefer and making the treatment-genotype combinations into individual groups. Then the contrasts are more obvious.

"The simplest manner is to not use a wald test, but rather an LRT" -- could you explain why?

Because treatment has multiple levels, so with a wald test one would need to handle the contrasts appropriately, whereas with an LRT one can simultaneously omit all the levels at once.


text:
Run cuffcompare in strand-agnostic mode
Is there a way to run Cufflinks' cuffcompare in a strand-agnostic mode?
I would like to do this because I have some RNA-seq datasets derived from an unstranded run, that should be compared to a reference transcriptome. I think that it makes no reason to have some transcripts in the experimental dataset classified with transfrag code tag x and s as compared to the reference if the input dataset was not derived from stranded protocol.
From manual:

x : Exonic overlap with reference on the opposite strand
s : An intron of the transfrag overlaps a reference intron on the opposite strand (likely due to read mapping errors)

Conceptually, does it make sense?
If so, and if the run by cuffcompare cannot be changed to a strand-agnostic manner, how should I interprete those transcripts classified as x? I was thinking of merging them to o and/or c .
Again, from manual:

c :   Contained
o :   Generic exonic overlap with a reference transcript


text:
Circos plot issue
I have taken out the data like chromosome number, start, and end coordinate then gene name but the circos need the data in a specific format as given in their tutorial circos
So far I have my data I have collected is as such
Symbol  Chr Start   End
RBM11   hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21;hs21 14216130;14216145;14216153;14216154;14216178;14219553;14219553;14219563;14219563;14221097;14221097;14221097;14221097;14221097;14224375;14224438;14224438;14224438;14224438;14226859;14226880;14226880;14226880;14226880 14216282;14216282;14216282;14216282;14216282;14219725;14219725;14219725;14219725;14221453;14221169;14221169;14221169;14221169;14224537;14224537;14224537;14224537;14224537;14227054;14228372;14227384;14228173;14228372
ELAVL4  hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1;hs1 50048014;50103910;50105690;50106292;50108933;50108966;50109626;50109706;50109816;50141772;50144953;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177130;50179544;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50200809;50200809;50200851;50200851;50200851;50200851;50200851    50048182;50104000;50105950;50106393;50109198;50109198;50109964;50109964;50109964;50141858;50145197;50145178;50145197;50145197;50145197;50145197;50145197;50145184;50145197;50145197;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50181300;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50197585;50197467;50197467;50197467;50197467;50197467;50197467;50197467;50203054;50201402;50202294;50203786;50201399;50203770;50201868
PTRF    hs17;hs17   42402452;42422627   42405388;42423517

So i have to make the the data as such
chr start end value [options]
...
hs3 196000000 197999999 71.0000
hs3 198000000 199999999 57.0000
hs4 0 1999999 28.0000
hs4 2000000 3999999 40.0000
hs4 4000000 5999999 59.0000

As of now i have multiple
RBM11   hs21;hs21;hs21;hs21;hs21;hs21

in for a single gene so how do i make it one ?or
RBM11   hs21

Any help or suggestion would be highly appreciated
I tried:
 sed -r ':a; s/\b([[:alnum:]]+)\b(.*)\b\1\b/\1\2/g; ta; s/(, )+/, /g; s/, *$//'

Updated answer so far
  library(splitstackshape)
df <- LEUKEMIA_CIRCOS_gene
Chr <-unique(cSplit(df, "Chr", sep = ';', direction = "long"))

this gives me this output
RBM11
hs21
14216130;14216145;14216153;14216154;14216178;14219553;14219553;14219563;14219563;14221097;14221097;14221097;14221097;14221097;14224375;14224438;14224438;14224438;14224438;14226859;14226880;14226880;14226880;14226880
14216282;14216282;14216282;14216282;14216282;14219725;14219725;14219725;14219725;14221453;14221169;14221169;14221169;14221169;14224537;14224537;14224537;14224537;14224537;14227054;14228372;14227384;14228173;14228372
2
ELAVL4
hs1
50048014;50103910;50105690;50106292;50108933;50108966;50109626;50109706;50109816;50141772;50144953;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50144957;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177089;50177130;50179544;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50193765;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50195561;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50197429;50200809;50200809;50200851;50200851;50200851;50200851;50200851
50048182;50104000;50105950;50106393;50109198;50109198;50109964;50109964;50109964;50141858;50145197;50145178;50145197;50145197;50145197;50145197;50145197;50145184;50145197;50145197;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50177192;50181300;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50193918;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50195786;50197585;50197467;50197467;50197467;50197467;50197467;50197467;50197467;50203054;50201402;50202294;50203786;50201399;50203770;50201868

my chromosome name are sorted out but the start and end remains because if I make concatenate the column I want to separate I get a really long list
    Start <- unique(cSplit(df,c("Chr","Start","End"),sep = ';', direction = "long"))

RBM11
hs21
14216130
14216282
2
RBM11
hs21
14216145
14216282
3
RBM11
hs21
14216153
14216282
4
RBM11
hs21
14216154
14216282
5
RBM11
hs21
14216178
14216282
6
RBM11
hs21
14219553
14219725
7
RBM11
hs21
14219563
14219725
8
RBM11
hs21
14221097
14221453

Thanks for providing a nice set of data to test with. However, could you tell what have you tried ans why it didn't satisfy it? We can then avoid suggesting the same route or improve your existing path.

I tried something like this from a solution in stack  sed -r ':a; s/\b([[:alnum:]]+)\b(.*)\b\1\b/\1\2/g; ta; s/(, )+/, /g; s/, *$//'  but i see trails of hs20;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;, after the hs20 as an example how to resolve that?

Well you seem to try to use a single regexp to solve a problem that needs to divide a line into several ones. I would count the ; and then pair each value of each column into a new line.

okay let me try to get around if i can do it myself ..

That's the spirit! Don't hesitate to come back if you have problems, or meanwhile perhaps I manage to answer it.

partially managed to do it through  a R library..i will update my question and you can see and suggest..

How are you reading this file? It seems like tab separated, but you seem to relay on cSplit to split by columns


text:
Database of position weight matrices for protein motifs?
I am trying to identify proteins that carry a consensus target sequence for a kinase. Usually if I am working with, for e.g. transcription factors binding to DNA I would use position weight matrices (PWMs) from one of a number of databases of PWMs to search the sequences with something like FIMO.
However, I don't know of any databases that contain protein motif PWMs. Specially I am looking for one for yeast Cdc28, but it seems like it would be a generally useful resource that I can't believe doesn't exist. For example, while the MEME suite was originally developed to identify protein motifs, and ouputs PWMs, all of its databases refer to DNA or RNA.

As the sequence of a protein usually is set on the DNA level, why should there be one specific for proteins? Given a DNA sequence you can find the amino acids, but not the other way around, so it is usually done with DNA (AFAIK)

Thats true, but usually DNA/RNA PWMs encode the binding sites of proteins on DNA/RNA. I am looking for the binding sites of proteins on other proteins. A DNA PWM is probably not suitable for this because, for example  Cdc28 binds the amino acid consensus S/T-P-X-K/R and there are many many ways of encoding this in DNA because each amino acid can be coded for by many codons. I doubt that a PWM built on the DNA sequences of the genes with this consensus would return a sensible PWM.

I think it easy enough to do it yourself with the sequences you want. You can use Biopython to build the PWM once you have the sequences (which could be downloaded from Uniprot).

My problem is I don't know the PWM - I only know the consensus sequence, not the frequency variation about that. That is why I was hopeing someone could point me to a database that might contain the full PWM. (e.g. The first amino acid of the target site can be Serine or Threonine, but I have no idea if one is favoured over the other)

And how has you built the consensus sequence? If it is consensus you have a MSA, so you can use it to calculate the PWM.

I just pulled it out of a paper. Actually, not true. My collaborator did or at least its something they just know and i assume they could look out the reference for if necessary...

PFam has HMM models for protein domains, this is almost certainly what you're after. You can download all PFam-A models from ftp://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz.
Then, you're looking for the "Pkinase", PF00069 family (I found this simply by searching the PFam website for "kinase"), so extract it from the file:
zcat Pfam-A.hmm.gz | awk -vRS='//' -vORS='//' '/PF00069/' > PF00069.pfam

Then, format it for HMMER format (you will need to install HMMER)
hmmpress PF00069.hmm

And now you can use it to search your sequences (fasta or flat file, probably others; see the documentation for details):
hmmscan PF00069.hmm prots.pep

I tried this on a fasta file with the sequences of CBK1_YEAST, CSK21_YEAST (kinases) and FLNA_HUMAN (non-kinase control) and got:
# hmmscan :: search sequence(s) against a profile database
# HMMER 3.1b2 (February 2015); http://hmmer.org/
# Copyright (C) 2015 Howard Hughes Medical Institute.
# Freely distributed under the GNU General Public License (GPLv3).
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# query sequence file:             prots.pep
# target HMM database:             PF00069.hmm
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

Query:       sp|P53894|CBK1_YEAST  [L=756]
Description: Serine/threonine-protein kinase CBK1 OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CBK1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    2.7e-64  203.2   9.0    3.5e-46  143.8   0.2    3.1  3  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 ?   -2.9   3.2      0.19      0.19     198     226 ..     222     243 ..     191     277 .. 0.57
   2 !  143.8   0.2   3.5e-46   3.5e-46       1     153 [.     352     505 ..     352     522 .. 0.95
   3 !   63.6   0.0   9.8e-22   9.8e-22     154     264 .]     568     672 ..     551     672 .. 0.83

  Alignments for each domain:
  == domain 1  score: -2.9 bits;  conditional E-value: 0.19
                           TTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
               Pkinase 198 sgekgkekvekeldqlekilkilgetkek 226
                           + +++++++++      + ++i ++++++
  sp|P53894|CBK1_YEAST 222 Q-QQQQQQQQQ------QHMQIQQQQQQQ 243
                           1.111111111......122222221222 PP

  == domain 2  score: 143.8 bits;  conditional E-value: 3.5e-46
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHH.HHHHHHHHHHHHHSSTTB--EEEEEEESSEEEEEEE--TTEBHHHH CS
               Pkinase   1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkke.kkvlrEikilkklkhpnivrlyevfeekdelylvleyveggeladl 88
                           +++++ +G+G+fG+V  + +k+tgki+A+K++ k+++ kk++  +v +E  +l+ +++p +v ly++f++ ++lyl++e+++gg+l+++
  sp|P53894|CBK1_YEAST 352 FHTVKVIGKGAFGEVRLVQKKDTGKIYAMKTLLKSEMYKKDQlAHVKAERDVLAGSDSPWVVSLYYSFQDAQYLYLIMEFLPGGDLMTM 440
                           67899********************************999999********************************************** PP

                           HHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEETTTEEEE-SGTTSEECCSSCS CS
               Pkinase  89 lsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLidekgelKitDFGlakelesssk 153
                           l r + ++e+ +++++++ + ++e +H+ g+iHrD+Kp+NiLid +g++K++DFGl++ ++++++
  sp|P53894|CBK1_YEAST 441 LIRWQLFTEDVTRFYMAECILAIETIHKLGFIHRDIKPDNILIDIRGHIKLSDFGLSTGFHKTHD 505
                           **999*****************************************************9999876 PP

  == domain 3  score: 63.6 bits;  conditional E-value: 9.8e-22
                           BCTCCSCGGGS-HHHHTTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH......TTS-HHHHHHHH CS
               Pkinase 154 ltsfvgtreYlAPEvlkeneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkeklpeaselseeakdllk 242
                           + s vgt++Y+APE+  +++y++++D+WslG i+ye l g ppf +e+ +e+++k+++          e+++++p + ++s ea dl++
  sp|P53894|CBK1_YEAST 568 AYSTVGTPDYIAPEIFLYQGYGQECDWWSLGAIMYECLIGWPPFCSETPQETYRKIMNF---------EQTLQFPDDIHISYEAEDLIR 647
                           6789*****************************************66665533333333.........557799*************** PP

                           HHT-SSGCCSTT....HHHHHTSGGG CS
               Pkinase 243 kllkkdpkkRlt....aeellqhpyl 264
                           +ll+    +Rl+    a+e+++hp++
  sp|P53894|CBK1_YEAST 648 RLLTHA-DQRLGrhggADEIKSHPFF 672
                           ***875.678877779********98 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (756 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.00
# Mc/sec: inf
//
Query:       sp|P15790|CSK21_YEAST  [L=372]
Description: Casein kinase II subunit alpha OS=Saccharomyces cerevisiae (strain ATCC 204508 / S288c) GN=CKA1 PE=1 SV=1
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------
    7.5e-62  195.1   7.2    3.2e-55  173.4   0.7    2.1  2  Pkinase   Protein kinase domain

Domain annotation for each model (and alignments):
>> Pkinase  Protein kinase domain
   #    score  bias  c-Evalue  i-Evalue hmmfrom  hmm to    alifrom  ali to    envfrom  env to     acc
 ---   ------ ----- --------- --------- ------- -------    ------- -------    ------- -------    ----
   1 !   25.5   1.1   4.3e-10   4.3e-10       1      60 [.      40      94 ..      40     115 .. 0.87
   2 !  173.4   0.7   3.2e-55   3.2e-55      56     264 .]     128     363 ..     101     363 .. 0.92

  Alignments for each domain:
  == domain 1  score: 25.5 bits;  conditional E-value: 4.3e-10
                           EEEEEEEEEESSEEEEEEEETTTTEEEEEEEEEHHCCCHHHHHHHHHHHHHHHHHSSTTB CS
                Pkinase  1 yekleklGeGsfGkVykaveketgkivAvKkikkekakkkkekkvlrEikilkklkhpni 60
                           ye ++k+G+G++++V+++v+ +++ ++++K++k  k+kk k     rEikil  l+ +++
  sp|P15790|CSK21_YEAST 40 YEIENKVGRGKYSEVFQGVKLDSKVKIVIKMLKPVKKKKIK-----REIKILTDLSNEKV 94
                           78999****************************99888855.....7******9998877 PP

  == domain 2  score: 173.4 bits;  conditional E-value: 3.2e-55
                            S.STTB--EEEEEEE..SSEEEEEEE--TTEBHHHHHHHHSS--HHHHHHHHHHHHHHHHHHHHTTEE-SS-SGGGEEEE.TTTEEEE CS
                Pkinase  56 k.hpnivrlyevfee..kdelylvleyveggeladllsrkgslseeeakkiakqilegleylHskgiiHrDlKpeNiLid.ekgelKi 139
                            + h ni++l++++++  +++  lv+eyv++ ++  l+    +l++ e+++++ ++l++l+y+Hs+gi+HrD+Kp+N++id ++++l++
  sp|P15790|CSK21_YEAST 128 NgHANIIHLFDIIKDpiSKTPALVFEYVDNVDFRILYP---KLTDLEIRFYMFELLKALDYCHSMGIMHRDVKPHNVMIDhKNKKLRL 212
                            35*************8878888*********9988887...5**************************************66677*** PP

                            -SGTTSEECCSSCSBCTCCSCGGGS-HHHH.TTSTCSHHHHHHHHHHHHHHHHHSS-STTTSSCHHHHHHHHHTHHHHHHHCTSHH.. CS
                Pkinase 140 tDFGlakelessskltsfvgtreYlAPEvl.keneyskkvDvWslGvilyelltgkppfsgekgkekvekeldqlekilkilgetkek 226
                            +D+Gla+ ++ + +++ +v++r +  PE+l ++++y+ + D+Ws G +l++++++++pf   +g++    ++dql+ki+k+lg+++++
  sp|P15790|CSK21_YEAST 213 IDWGLAEFYHVNMEYNVRVASRFFKGPELLvDYRMYDYSLDLWSFGTMLASMIFKREPFF--HGTS----NTDQLVKIVKVLGTSDFE 294
                            ******************************999***************************..3443....789**********77766 PP

                            ....TTS...............................-HHHHHHHHHHT-SSGCCSTTHHHHHTSGGG CS
                Pkinase 227 lpeasel...............................seeakdllkkllkkdpkkRltaeellqhpyl 264
                                +                                 ++e++dl+++ll++d+++Rlta+e++ hp++
  sp|P15790|CSK21_YEAST 295 KYLLKYEitlprefydmdqyirkpwhrfindgnkhlsgNDEIIDLIDNLLRYDHQERLTAKEAMGHPWF 363
                            5554444588999999****************************************************9 PP

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (372 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         1  (1); expected 0.0 (0.02)
Passed bias filter:                        1  (1); expected 0.0 (0.02)
Passed Vit filter:                         1  (1); expected 0.0 (0.001)
Passed Fwd filter:                         1  (1); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               1  [number of targets reported over threshold]
# CPU time: 0.01u 0.00s 00:00:00.01 Elapsed: 00:00:00.02
# Mc/sec: 4.91
//
Query:       sp|P21333|FLNA_HUMAN  [L=2647]
Description: Filamin-A OS=Homo sapiens GN=FLNA PE=1 SV=4
Scores for complete sequence (score includes all domains):
   --- full sequence ---   --- best 1 domain ---    -#dom-
    E-value  score  bias    E-value  score  bias    exp  N  Model    Description
    ------- ------ -----    ------- ------ -----   ---- --  -------- -----------

   [No hits detected that satisfy reporting thresholds]

Domain annotation for each model (and alignments):

   [No targets detected that satisfy reporting thresholds]

Internal pipeline statistics summary:
-------------------------------------
Query sequence(s):                         1  (2647 residues searched)
Target model(s):                           1  (264 nodes)
Passed MSV filter:                         0  (0); expected 0.0 (0.02)
Passed bias filter:                        0  (0); expected 0.0 (0.02)
Passed Vit filter:                         0  (0); expected 0.0 (0.001)
Passed Fwd filter:                         0  (0); expected 0.0 (1e-05)
Initial search space (Z):                  1  [actual number of targets]
Domain search space  (domZ):               0  [number of targets reported over threshold]
# CPU time: 0.00u 0.00s 00:00:00.00 Elapsed: 00:00:00.00
# Mc/sec: inf
//

HMMER is a great tool and can also analyse multiple sequence alignments to build the hmms if that's what you need. I strongly urge you to reads its documentation which is both informative and a surpisingly fun read (it includes such phrases as "Their format is “proprietary”, which is an open source term of art that means both “I haven’t found time
to document them yet” and “I still might decide to change them arbitrarily without telling you”").

I was sure that with pfam and hmmer it could be possible to find a specific protein motif. But given that the OP has a consensus it might be hard to use it to search in pfam.

@Llopis pfam has HMM models which are essentially PWMs. And this isn't a way of searching PFam, it's a way of using PFam to search protein sequences. The OP says: "I am trying to identify proteins that carry a consensus target sequence for a kinase.", so I assume Ian has a list of proteins and needs to search through them for a kinase domain. I may be wrong, but that's what I got from the question.

I need to review my notes on HMM. Thanks for the answer!

Thanks, but its the pwm for the kinase targets rather than the kinase domain i was looking for.


text:
What is the best method to estimate a phylogenetic tree from a large dataset of >1000 loci and >100 species
I have a large phylogenomic alignment of >1000 loci (each locus is ~1000bp), and >100 species. I have relatively little missing data (<10%).
I want to estimate a maximum-likelihood phylogenetic tree from this data, with measures of statistical support on each node.
There are many phylogenetics programs that claim to be able to analyse datasets like this (e.g. RAxML, ExaML, IQtree, FastTree, PhyML etc). Given that I have access to a big server (512GB RAM, 56 cores), what are the pros and cons of each program.
Which programme likely to give me the most accurate estimate of the ML tree for a dataset of this size within a reasonable timeframe?

You want us to review all the phylogenetic programs or just the ones you listed: RAxML, ExaML, IQtree, FastTree, PhyML (Which still I think is quite broad). How do you measure the accuracy of the ML tree?

I'm interested in any and all opinions, evidence, and links to comparisons between any software capable of estimating ML trees from large datasets like this. There would be many ways of measuring accuracy, including: (i) evidence from simulation; (ii) comparisons of Likelihood scores of estimated trees in a common framework.

RAxML is the goto program. The problem with FastTree is that its accuracy is limited to the datasets used, i.e. yours might be an outlier.

This paper claims that FastTree is almost as accurate as RAxML, while being much faster. You just have to be careful, however, that the support values output by FastTree are not bootstrap values, they are based on the Shimodaira-Hasegawa test. (Also, see this comment for the case you have very short branch lengths).   [update: However, according to the recent comparison paper mentioned below FastTree performed quite poorly in comparison to RAxML or IQ-tree.]
From what I understand, you should use ExaML only if your data is too large to be handled by RAxML in a single node. ExaML should perform like RAxML but with some parallelization overhead. For all effects I treat them as the same. I don't know of relevant advantages of phyML over RAxML (for me, it's easier to use but I am very used to phyML).
I am not familiar with IQ-tree, but its authors claim that even given the same time as RAxML or phyML, IQ-tree already finds better likelihoods more often than not (although by default it takes a bit longer to converge). A recent comparison between all these programs favoured IQ-TREE for both single-gene and concatenation analyses (with RAxML very close). It may also estimate branch support through a SH-like test only, but I'm not sure.   [update: IQ-tree offers 3 measures of support, standard bootstrap, aLRT, and ultrafast bootstrap. See OP's comment below for details.]
However, since you have few missing data, you might also want to try a single-locus tree inference followed by gene tree clustering (using treescape or treeCL) to see how spread your data is, or to see the effect of removal of outliers, or to use ideas similar to statistical binning.

In the recent comparison you mention, I note that RAxML only performed better than IQtree when they did 10 independent searches per replicate. So one might just favour IQtree here too, noting that it's sensible to do a bunch of independent searches. Also FastTree performed very poorly in all comparisons of that ms.

Measures of support in IQtree: (i) standard bootstrap; (ii) aLRT (approximate likelihood ratio test, which is something like asking whether a given branch length is significantly >0); (iii) Ultrafast bootstrap (not the same interpretation as a standard bootstrap, more like a posterior probability if I've understood it right).

Note that likelihood comparisons between program outcomes also depend on the choice of model and data partitioning. I haven't looked at the details, but some programs may implement models that are not available in others. This can be a criterion of choice.

Thank you for the comments @roblanf, I took the liberty of updating my answer with this info.

note also that there is a newly released re-write of RAxML available here. It combines the best parts of RAxML and ExaML, while being faster in general. However not all current RAxML features have been implemented yet.


text:
How do I carry out an ancestry/admixture test on a single VCF file?
This is a question from /u/beneficii9 on reddit. The original post can be found here.

Through the Personal Genome Project, I have had my whole genome sequenced by Veritas, and have it in the form of a single VCF file for the whole genome and one BAS file for each chromosome. The reference genome associated with the VCF file is hg19. It has been helpful in health data; for example, I discovered I'm homozygous for the non-functional variant CYP-2D6 gene (rs3892097), which can render several common medications useless, and helps explain why some medicines didn't really work for me. My doctor has found this information very helpful.
Unfortunately, I can't find any way of looking at admixture or ancestry. I've tried setting everything up using a combination of VCFTools, Plink1.9, and ADMIXTURE, but I can't get it to work. I think for ADMIXTURE you have to have a bunch of genomes sorted by geographical origin to compare your genome against, but I'm not sure how to do that, and what's online isn't very clear to me. So scratch that one off.
I've tried converting the file to 23andme format (and at this /u/psychosomaticism has been very helpful). I did that (though it seems there were problems because of the way the VCF file was set up). But the websites that take the data want you to point them to your 23andme account, and that doesn't really work if you only have the file. 23andme doesn't provide for people who had their whole genomes sequenced. They want you to give them a saliva sample like everyone else.
So, what can I do?

A modified implementation of Vivek's answer.
peddy is a Python package that samples an input .vcf at ~25000 sites and projects onto a principal component space built on 2504 thousand genome samples. The author has extensive documentation of the tool's features and a link to the preprint.
I downloaded the .vcf and .vcf.tbi for the NA12878 sample from Genome in a Bottle's ftp here. Then, created a custom .ped file, NA12878.ped, with the contents:
NA12878 HG001   0       0       2       0
At the command line:
$ peddy --plot --prefix myvcf HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz NA12878.ped

The output files all have the prefix myvcf., here's myvcf.pca_check.png

I think it might be worth adding: peddy requires coloredlogs and cyvcf2 the later of which will required you to have gcc&g++ and the headers (e.g. -dev versions in ubuntu) of libbz2, liblzma and curl. This knowledge would have saved a precious 15 precious moments of my life :)

Also, now that it's installed I seem to get the error: AssertionError: error loading tabix index for b'HG001_GRCh37_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_highconf_PGandRTGphasetransfer.vcf.gz' (Also based on my understanding shouldn't the 2 be 0 or left to the user to set based on his gender ?)

A suggestion from the UK Biobank QC methods white paper:

Create a set of SNPs common to your VCF and the 1000 genomes phase 3 call set.
Perform PCA of the 1000 genomes samples using eigenstrat smartpca. You might have to convert to binary plink format.
Project your genotypes onto this pre-computed PCA space and visualize using ggplot to see which cluster you fall in.

Ancestry testing is a tricky subject. I spent a good chunk of my PhD project on questions similar to this, and haven't really found a good answer for how to detect unmodelled ancestry.
The general idea of how I tried to approach ancestry determination was to create a model set containing well-defined groups of individuals with a known, specific ancestral background. A query individual (or individuals) were then added to this group, and an ancestry-estimating program was run (specifically Structure) to work out what proportion of the unknown individual was attributable to each of the known groups. I would expect that most genetic ancestry tests follow a similar approach to this, although perhaps with a little bit less care about the definition of the model population groups.
Problems arise when individuals are incorrectly assigned to a particular group, when groups are present in the model set that represent a larger proportion of individuals than most other groups, when groups are closely-related to other groups, and when an ancestral history is present in a test individual that doesn't match any of the model groups. And all that assumes that the marker set used for ancestry determination is perfect: no bias towards any particular group, and no systematic genotyping error.
This doesn't mean that ancestry testing won't work, but it's a good idea to take the results with a large grain of salt. There's a good example of a media personality in New Zealand who was told that she had a very high probability of being 100% Māori, despite having a good knowledge of her own family history that indicated one European ancestor a few generations back on both sides of her family.

Adding another answer that, while much more involved, offers flexibility, walks through minimal steps, and provides visualizations.
I created a Dash app and repository that explores this type of analysis using published ancestry informative SNPs and 1000 Genomes Project Data.

Identify published ancestry informative SNPs (AISNPs).

Kidd et al.
Kosoy et al.

Download 1000 Genomes Project genotype bcf data.
Limit genotype data from (2) to AISNPs loci from (1).
One-hot encode the genotypes.
Perform dimensionality reduction (PCA, t-SNE, or UMAP)
Plot components.

The original post seemed to want to project their genome against a reference sample. The plot_walkthrough.ipynb in the tgviz repo shows how this can be achieved.
Output from the notebook on my own data:

Related: ethnicity check either from bam or vcf files
The missing piece of your puzzle is an ancestry-labeled reference dataset.  1000 Genomes phase 3 usually works well enough here, and it now has variant callsets for both GRCh37 and GRCh38.


text:
Finding data to validate results
My group has a complex datset and I would like to validate the results and check the methods I use in a similar dataset (better if it has been already studied).
Characteristics of my dataset:

It is related to the inflammatory bowel disease (IBD)
More specifically to Crohn's disease.

It has RNA-seq
I have explored the GEO and there are few projects (10) with RNA-seq on the tissue of my disease, and I don't know how to find the 16S-seq of those projects if it exists (see point below).

It has 16S-seq (and not only from the stools)
The ENA seems to have 16S sequencing but most of them are from stools not from biopsies and still, I don't know how to find the RNA-seq of those projects if it exists (see point above)

It has both RNA-seq and 16S-seq from the same individuals at the same time point.
The NIH Human Microbiome Project for my disease has only one individual with that both 16S-seq and RNa-seq of the same patient at the same time.

How can I find a dataset similar to mine?

The snarky answer is that it's better to find a matching dataset at the experimental design phase (i.e. before these studies are carried out), so that the bioinformatics component on the other side of things ends up being a lot easier and more robust. Although it can be quite time consuming to do a database search at the start, it's even harder to find matching projects if the legwork hasn't been done during the design phase. Researchers should be thinking about how hypotheses will be validated before the experiments are carried out.
Every additional characteristic that is chosen for comparison adds complications with regards to matching data, and reduces the applicability of results when there is a mismatch. You would have to be somewhat lucky to be able to find a IBD-Crohn's dataset with paired RNA-seq data and 16s-seq (presumably short-read; is 454 good enough?) from a specific non-stool tissue.
Apart from the places you've already looked, there's dbGaP and SRA. Given the number of characteristics that have been mentioned, I expect you'll have to make some compromises when comparing other datasets.

Thanks for the tip, but is hard to find projects at the experimental design phase. We have an easier dataset in our lab, but I would prefer to confirm it with other datasets. I had forgotten about these databases, thanks.


text:
arm64 installation of AMRFinderPlus
I was trying to install NCBI AMRFinderPLUS on my Nvidia Jetson using standard Bioconda based installtion. However, all the time it is giving error that the package is not available and I thought I am making some typo or channel is not proper.
(amrfinder) jay@jetmimd-desktop:~$ conda install -y -c conda-forge -c bioconda ncbi-amrfinderplus
Channels:
 - conda-forge
 - bioconda
 - defaults
Platform: linux-aarch64
Collecting package metadata (repodata.json): done
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - ncbi-amrfinderplus

Current channels:

  - https://conda.anaconda.org/conda-forge
  - https://conda.anaconda.org/bioconda
  - defaults

To search for alternate channels that may provide the conda package you're
looking for, navigate to

    https://anaconda.org

and use the search bar at the top of the page.

Then I remembered about the different architecture. So, now here I am. Is there any way to install this on arm64 architecture?
Thank you for the help in advance!!!!

Try the following
conda install bioconda::ncbi-amrfinderplus

The above does work I tested it (below). There could be an issue with you version of conda (24.1.2 OSX).
Also has the package has been ported for GPU(?) ... I assume it is, otherwise it's obviously not going to work and conda will not install it
However, if it ported for GPU and still fails use mamba (I don't use mamba BTW), because it is very slow and that could be where the error is (unusually slow).
My connection initially was unsuccessful but resolved via

Retrying with flexible solve.

To demonstrate I got it working here's the reply ...
The following NEW packages will be INSTALLED:

  blast              bioconda/osx-64::blast-2.12.0-pl5262h78c34c6_0
  c-ares             pkgs/main/osx-64::c-ares-1.19.1-h6c40b1e_0
  curl               pkgs/main/osx-64::curl-7.88.1-hdb2fb19_2
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  hmmer              bioconda/osx-64::hmmer-3.3.2-hb19a6fa_3
  libcurl            pkgs/main/osx-64::libcurl-7.88.1-ha585b31_2
  libev              pkgs/main/osx-64::libev-4.33-h9ed2024_1
  libidn2            pkgs/main/osx-64::libidn2-2.3.4-h6c40b1e_0
  libnghttp2         pkgs/main/osx-64::libnghttp2-1.52.0-h1c88b7d_1
  libssh2            pkgs/main/osx-64::libssh2-1.10.0-hdb2fb19_2
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  ncbi-amrfinderplus bioconda/osx-64::ncbi-amrfinderplus-3.11.14-hfee4e81_0
  perl               pkgs/main/osx-64::perl-5.26.2-h4e221da_0
  perl-archive-tar   bioconda/osx-64::perl-archive-tar-2.32-pl526_0
  perl-carp          bioconda/osx-64::perl-carp-1.38-pl526_3
  perl-common-sense  bioconda/osx-64::perl-common-sense-3.74-pl526_2
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.087-pl526h6de7cb9_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.087-pl526h770b8ee_0
  perl-exporter      bioconda/osx-64::perl-exporter-5.72-pl526_1
  perl-exporter-tiny bioconda/osx-64::perl-exporter-tiny-1.002001-pl526_0
  perl-extutils-mak~ bioconda/osx-64::perl-extutils-makemaker-7.36-pl526_1
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.087-pl526h6de7cb9_0
  perl-io-zlib       bioconda/osx-64::perl-io-zlib-1.10-pl526_2
  perl-json          bioconda/osx-64::perl-json-4.02-pl526_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl526h04f5b5a_3
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-0.428-pl526_1
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.428-pl526_0
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl526h1de35cc_1
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.52-pl526h01d97ff_0
  perl-types-serial~ bioconda/osx-64::perl-types-serialiser-1.0-pl526_2
  perl-xsloader      bioconda/osx-64::perl-xsloader-0.24-pl526_0
  wget               pkgs/main/osx-64::wget-1.21.4-ha585b31_1

There is a blast gpu and I'd assume the other library/packages have gpu versions. There's a lot of Perl code there BTW, so it's an old package, before NVIDIA was big perhaps?

If there is a GPU compliant channel on conda (nanoporetch) - that is impressive: especially given the amount of Perl code. Porting for GPU from CPU isn't trivial.

Few versions behind

The key thing is whether the upgrades are fixing bugs or refactoring/adding functionality. If it's the latter then I would go  for it because the GPU would compensate for the slower runtime (obviously given you don't need the extra functionality).

Hey, Thank you very much! I had tried that yesterday  (among the many other things I could try/knew how to try) after posting here. But that also generated extact same error as posted above. I think(it may be basic and I am still not sure) the packages which work for OSX64 may not work for aarch64/arm64 and that is why I am not able to use this streamline way of installation. Although I also came across, separate channel on conda nanoporetch who have stated many packages for aarch64, just one or two versions behind. So, I am still in dilemna what to do!  Use the older versions or ...

infact I am not able to install the packages with noarch tag. Now, I am not sure weather this is a trivial information to know how to work around these packages or this is something which will take time in developing.

@prekij if there's no available package, you will probably have to compile it yourself for the architecture :/


text:
Identification of metabolites through spectral library search
I am very new to single-cell metabolomics and today I want to jump in this area. I have been reading many review papers related to identification of metabolites at single-cell resolution from spectral data through spectral library search and I have a question that cannot find out a solution.
The common workflow to identify mebolites I have discovered so far is that starting from, let say, 1000 live cells, we extract the content of each cell, put them to, let say, GC-TOF-MS technology to have spectral graphs (x-axis presents m/z values and y-axis presents intensities), and finally compare the spectral graph of each cell against a certain spectral library to infer metabolites.
The thing is we know that a single cell has multiple metabolites regulating cellular functions of that cell. However, the data I see show that each cell has only one spectrum (the figure below) => we can only find out 1 metabolite for one cell??? Is that true? Or something I am missing?

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


text:
(Complete) Stuctural Models of Antibodies: Combining Fc and Fab from pdb to Build Hinge Region? Or what?
I'm a newcomer to the field of antibodies, and I'm seeking guidance on obtaining 3D structures for my molecular simulations. After browsing the PDB databank, I noticed that there are Fab and Fc structures available. My specific interest lies in investigating the structural dynamics in the entire antibody (AB) and on the hinge region.
I have a few questions:

Is it common practice to attempt building the entire antibody structure and hinge regions from PDB files for fragments, or are there recommended alternatives?

Are there comprehensive PDB structures available for different IgG and therapeutic antibodies that include information on the whole structure and hinge regions? Especially for Bevacizumab and Rituximab?

Does it make sense to select single segments, e.g., CH1, VH, VL, for analysis? Would one expect these structures to be reasonably stable?

Additionally, I'd appreciate any insights into common pitfalls and any recommendations you might have for someone new to this field.
Thank you in advance for your time and assistance.
Best regards,
Albert

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

Hi, if you've got multiple questions, it's often better to split them up into separate question posts. However, exploratory posts like this are almost always going to produce better outcomes on discussion websites like reddit

wait for cryomicroscopist to deliver https://www.nature.com/articles/s41467-022-34090-2  Here, we use cryo-EM to image the structure of the human full-length IgM pentamer, r


text:
Suggestion for way foward with upregulated genes
I am working with single cell RNA data- control and treated. I analyzed and got a list of upregulated genes. Upon doing pathway analysis all the genes were of apoptotic pathway. But these genes can signify different meaning to my interpretation, because, these genes were upregulated even in diseased condition. let me take an example, NEAT1 is one of the gene that is seen upregulated in the cancer condition and is associated with cancer progression. the treated sample which is a therapy also has neat1 upregulated compared to control which was untreated. Just because pathway analysis showed apoptosis, I cannot conclude that neat1 gene is aiding in killing the cancer cell due to the therapy right/ It can also be like due to continued upregulation its delaying cancer cell death due to therpay.
How to confirm computationally the above two statements and find the upregulated genes actual role in the therapy?
to be noted: there is no prior data or pattern available for the therapy based gene expression yet.


text:
Kallisto RNA-seq errors because of incompatible indices when running <kallisto quant> using downloaded mouse index
I'm trying to run RNA-seq on single-end mouse RNA data (50bp reads) but have very limited RAM (<4 gb) - I don't think I really have a way to work around this memory restriction, and wanted to try Kallisto's pseudoalignment (I conda installed it) as indexing and/or aligning with other tools like STAR and HISAT2 have resulted in the processes being killed. Building the index with kallisto index still results in the process being killed, so I downloaded the index files from https://github.com/pachterlab/kallisto-transcriptome-indices/releases/tag/v1 instead. I ran kallisto quant -i index.idx -o output --single -l 180 -s 20 sample1.fastq (I don't entirely know what the mean and sd arguments should be for my data), but this resulted in the following:
[quant] fragment length distribution is truncated gaussian with mean = 180, sd = 20
Error: incompatible indices. Found version 13, expected version 10
Rerun with index to regenerate

I've seen that other people have come across this error before when quantifying with kallisto, but it seems that they are building the index themselves while I was trying to use the pre-generated index from the unpacked mouse_index_standard.tar.xz file downloaded from the link above.
How can I get around this error? And in general, are there any suggestions on how to get raw/normalized count data from RNA-seq data with such memory limitations? I'm rather new to this.

If your finances at all allow, get more RAM. You're not soing much with such limited resources in the long run.

It means that the version of the tool used for indexing and your current version are different and not compatible. Either downgrade your version to the version used for indexing, or build a new index with your current version. As by the GitHub repo, these indices that you link (supports kallisto versions >= 0.50.1), so use this version if you want to use the index.
If that does not work then I would open an issue at their GitHub.


text:
How to download Sequencing data on Windows using SRA toolkit?
I have downloaded and installed the SRA toolkit, but there seems to be no online material on how to download SRR files using the SRA toolkit on Windows.
All online materials are Ubuntu or Linux based.
Is there some obvious resource that I am missing?

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

I found a solution to this problem.
After having installed SRA toolkit:
prefetch -v SRR9....1
this helps download the SRA file.
Then to extract it as FASTAQ
fasterq-dump SRR9....1

That just postpones the problem. Now you need to align these files, tell me a single popular Widows aligner.

@ATpoint clustal? Please try and make your comments helpful instead of just snide. Of course you cannot seriously do bioinformatics if you just use Windows, but people who only dabble might be able to get away with it, and the specific example of alignment tools is a bad one.


text:
Prediction of DNA binding motif from sequence/structure of TF/regulator protein?
Do you know if it is actually possible to predict with good accuracy a DNA binding motif from only the sequence/structure of the DNA-binding protein?

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

This question looks like it's too general to be useful for asking or answering; please frame your question around a specific problem that you're having, especially including detail about the nature of the DNA binding motif that you're interested in.

Find papers....

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6153143/


text:
Best Library for 3D Protein Structure Visualization and Annotation
I have a background in Computer Science and have recently gotten very interested in the field of Protein Engineering. I am trying to build a web app which displays the 3D structure of a protein when I upload a PDB file. I also want to annotate subsequences with unique labels.
As I am very new to the field I am confused about which library to use. Can I use BioPython for everything? Or should I use any other library for 3D structure visualization like 3Dmol.js, NGL Viewer etc?
Any comments/suggestions/feedback would be really helpful guys as I am a beginner in the field of BioML.

BioPython can't do structure visualization. NGL and LiteMol are obsoleted by Mol* – from the same developers. There is venerable J(S)mol and its wiki lists a few other WebGL viewers. The newest viewer is Moorhen, but it's more for structural biologists.

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


text:
Find up/downregulated gene lists across scRNA-seq samples
I have two separate scRNA-seq datasets that are different batches. I want to get genelists of up/downregulation of one cluster in one dataset compared to another cluster in the other dataset.
I already did integration using scanorama on Scanpy and the anchor method on Seurat (I can use either tool) but the Satija Lab clearly says on github that you can't do differential expression analysis on batch-corrected data.
Can I just append anndata/seurat objects from different batches and do DE without batch correction? Or is there something I'm missing?


text:
Seeking Advice on Acquiring Clean Experimental data comprises of proteins for FDR Estimation
I am a bioinformatics student doing my research on top-down mass spectrometry-based proteomics and seeking some guidance on acquiring very clean experimental data with known true positives and false positives. I want to assess my proposed FDR estimation method under controlled conditions and would appreciate your insights on obtaining high-quality experimental data.
My question is

Any tips on acquiring clean experimental data with known true positives?
Any suggestion on sample types or protocols you recommend for ensuring a reliable dataset for FDR estimation?
How can one introduce controlled errors or variations, like false positives, while maintaining the overall data quality? If anyone has any idea or experience or any paper suggestion in which they used a similar approach please share thanks.

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


text:
Where can you download functional domain labels of the human genome?
I'm visualizing where mutations occur in a gene, and would like to add functional domains to the plots.
I have some trouble locating sources to download data on functional domains.
The dataset should contain:

Functional domain name
Start position in gene
End position in gene
Gene symbol name

Any tips kindly appreciated!

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

I would suggest looking into GFF3 annotations of the human genome.
If what you want is protein domain annotations, you may be more interested in the UniProt database which does seem to include genome annotation tracks for different protein functional regions.
I am unclear if uniprot directly has the gene names but this would be trivial to generate with a bedtools intersect operation.


text:
What is a common number of cells per cell type in single cell RNA Seq?
Good afternoon,
I just started working with Single cell RNA Seq and I am trying to understand what is an average / possible number of cells per sample.
I am looking at a dataset of 7 patients with kidney tumour samples (please find attached). Would these number of cells be possible? I basically want to double check whether I didnt do any mistakes in my pre-processing.
Thank you

It depends on many factors, both biological and technical, and it's impossible to answer in general.
The starting question is which technology was used, and how many cells were used in input. For example, in a SMART-Seq experiment with one cell per well (in a multi-well plate), you should know how many cells you have (one per well). For a 10x experiment, the person running the experiment should have measured the cell concentration in the suspension before loading the chip, and targeted a particular number of cells (typically 10,000). In a combinatorial barcoding approach (e.g. Parse), you'd also have a measurement beforehand so as to load a particular number of cells per well. So, you should already have an idea of the total number of cells.
Then, the first step of the analysis is typically some QC. In barcode-based approaches (including 10x), you'd make a knee-plot to estimate the actual number of sequenced cells in your sample.
In the clustering and annotation steps, you'd start distinguishing and annotating groups of cells, which hopefully do correspond to cell types. The algorithms are not perfect, so the clusters don't necessarily match cell types perfectly. For example, a rare cell type will easily be lumped together with another, similar, cell type. For a well-represented cell type, the algorithms have more chance to capture subtypes or cell states.
And of course there is the biology, and your knowledge of it. How many astrocytes are there in the kidney? I assume any number above 0 is worrisome. I'd expect the cell types in the data to roughly reflect the cell types in the tissue (but not exactly, there can be all kind of bias).
There is a question of depth: if you have few cells in a cluster, and few reads per cell, the gene expression data in that cluster is not very reliable. The actual number depends on the technology and read depth. I would say for 10x I would prefer at least 100 cells, and would not trust a cluster with <10 cells.
So for your question, you need to use your knowledge of the kidney and of the experimental design to decide. You can also compare to published scRNA-Seq experiments of the same tissue. You can also compare between patients to see if there is something troubling (for example, why does patient 1 have epithelial cells but no adipocytes, and opposite for patient 7? It might not be a problem, but only with context can you tell). And importantly, you do need to check that the annotation does match known gene expression patterns.

Dear Alexlok, Thank you so much for your detailed and very helpful answer, I will save it in my notes :). I have read a couple of times now the sentence "use your knowledge of biology" e.g. the found cell clusters. I definitely need to increase my cell biology knowlege here I think.... And yes also consider the experimental design more.