text:
Why is bulk RNA sequencing reflecting AVERAGE expression but not TOTAL expression of all cells?
When I am reading papers that compares bulk RNA sequencing and single-cell RNA sequencing, we often see papers describe bulk RNA seq measures the average cell expression.
For example, in this paper Single-Cell RNA-Seq Technologies and Related Computational Data Analysis

bulk RNA-seq mainly reflects the averaged gene expression across thousands of cells

And this paper A practical guide to single-cell RNA-sequencing for biomedical research and clinical applications

the averaging that occurs in pooling large numbers of cells does not allow detailed assessment of the fundamental biological unit—the cell—or the individual nuclei that package the genome

And this paper Differential expression analyses for single-cell RNA-Seq: old questions on new data

Standard RNA-Seq experiments need millions of cells for sequencing [2,3], and therefore can only get averaged measurements of gene expressions of the cells sequenced.

My question is, when we perform the bulk RNA seq and calculates parameters like TPM, we don't actually divide the total number of mRNA transcripts with the number of cells (I think sometimes we don't even have an accurate number on the number of cells that we used in bulk RNA seq), we simply normalise it with the transcript length and count.
Therefore, if we DID NOT divide the number of transcripts with the total cell number, how are we measuring the average of expression, but not the total expression?
Or in other words, what is the "average expression" referring to? Average to what?

What these statements refer to is that bulk RNA-seq will overlook the different transcriptional cell types, hence the value of single cell RNA-seq in unraveling heterogeneous cell populations/tissues.
I think "total" is also avoided because total RNA has another meaning related to capturing all the RNA transcripts in a cell (e.g. including non polyA transcripts), which is still not the case for most RNA-seq experiments.
Also note that RNA-seq experiments generally (without spikins) are only relative and not giving absolute quantifications of the amount of RNA in a cell or groups of cells.

Thank you for your answer! I understand what you said, but I'm asking about total expression, not total RNA. You said "bulk RNA-seq will overlook the different transcriptional cell types", yes I agree, but then it also means that it is interpreting the total transcriptional actives within the mixture of cells. I'm still not sure why the word "average" is used in these publications. Or in other words, what is the "average expression" referring to? Average to what?

@benson23 As I said, these numbers are generally relative and not absolute so these statements are not implying an actual mathematical calculation of mean expression (or indeed summed total expression). 'Average' also means 'typical' and that is the essence here, but this is semantics

Thank you for explaining again! I agree that these numbers would be relative, however, I do not quite agree with you that the use of the word "averaged" is semantic. Based on the first paper I quoted, it said "averaged gene expression across thousands of cells", in the second paper "the averaging that occurs in pooling large numbers of cells" and in the third paper "only get averaged measurements of gene expressions of the cells sequenced", which to me, certainly implies some mathematical operations...

Moreover, in academia, the words used should be precise. If the meaning here really refers to "typical", I believe these authors would not use such an ambiguous word of "average". If it really doesn't have a mathematical meaning, the phrase "relative expression" or as you said, "typical expression" would be much better. Also, if you try to replace the word "averaged" with "typical" in the examples that I quoted, it'd be kind of weird.

I believe using the average is fair for this purpose. The moment you break thousands or millions of cells you are already averaging relative mRNA content; you are making a soup of it. Keep in mind that when you sequence RNA (RNA-seq), you are not sequencing every single molecule in the sample. In other words, you do not know how many molecules of each mRNA you have.

What you are doing is comparing relative level of certain mRNAs between two samples in a giving sequencing depth. It's like comparing salinity iin two different seas: Dead Sea has salinity of 337 g/kg and Black Sea has salinity of 20 g/kg. You can consider the salinity a single mRNA. You do not need to know the exact size of Dead Sea to say that salt concentration of Dead Sea is ~16-fold higher than Black Sea.


text:
How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
I would like to download a number of SARS-CoV-2 fastas from NCBI/Genbank using a web endpoint. I looked on their website but it's confusing.
I have the accession numbers, but can't figure out how to use the NCBI/Genbank Entrez API.
This is how I would do it using cov-spectrum: https://lapis.cov-spectrum.org/open/v1/sample/fasta?genbankAccession=MZ039713,MZ039682,OV522182
Is there something similar for NCBI/Genbank?
There must be a way, I just can't figure it out. The whole entrez thing is overwhelming to me, all I want is a simple API, like covSpectrum offers, for example. I scanned this but didn't help: https://www.ncbi.nlm.nih.gov/sars-cov-2/
E.g. how could I download all these SARS-CoV-2 nucleotide sequences with accession numbers: MZ039713,MZ039682 and OV522182 using a URL?
Some command line tools are advertised, but I don't want to have to use a CLI tool. Just a simple URL, does that not exist? CLI tools: https://www.ncbi.nlm.nih.gov/datasets/docs/v1/quickstarts/command-line-tools/
I also don't want to do any clicking through a web interface, just a URL so that other web apps can download it behind the scenes.

Could you please let us know which ones?

@M__ what do you mean by "which ones"? For example fastas for these accession IDs: MZ039713, MZ039682 and OV522182 like displayed here: https://www.ncbi.nlm.nih.gov/nuccore/OV522182

Thanks I'll put a response together. Can't do it immediately, but will do today

You can get the data from the API of the European Nucleotide Archive (ENA) which is synched with NCBI/GenBank:
https://www.ebi.ac.uk/ena/browser/api/fasta/MZ039713,MZ039682,OV522182

You can try the NCBI Datasets API. We have an endpoint that would allow you to download virus sequences by accession  The only thing to be aware (based on your list of accessions) is that our API needs the accession version to return results. In the case of your examples, they would need to be MZ039713.1, MZ039682.1 and OV522182.1.
I hope this helps!

Thanks! Can you provide an example query? The API documentation is a bit overwhelming to me. Also, am I correct in assuming you need a private API key?

You don't need a private API key.
Using the accessions you provided as an example, you can create a json file (accessions.json):
{"accessions":["MZ039713.1", "MZ039682.1", "OV522182.1"]}
Then, you can use the following command to download everything using the API:
cat accessions.json | curl -H "Content-Type: application/json" -X POST --data-binary @- "https://api.ncbi.nlm.nih.gov/datasets/v1/virus/genome/download" --output file.zip
Let me know if you run into any issues.

In the efetch.fcgi URL, you can set rettype=fasta to download the nucleotide sequence or rettype=fasta_cds_aa to download protein coding sequences:
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta&id=NC_045512.2'|head -n3
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC
$ curl 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&rettype=fasta_cds_aa&id=NC_045512.2'|head -n3
>lcl|NC_045512.2_prot_YP_009724389.1_1 [gene=ORF1ab] [locus_tag=GU280_gp01] [db_xref=GeneID:43740578] [protein=ORF1ab polyprotein] [exception=ribosomal slippage] [protein_id=YP_009724389.1] [location=join(266..13468,13468..21555)] [gbkey=CDS]
MESLVPGFNEKTHVQLSLPVLQVRDVLVRGFGDSVEEVLSEARQHLKDGTCGLVEVEKGVLPQLEQPYVF
IKRSDARTAPHGHVMVELVAELEGIQYGRSGETLGVLVPHVGEIPVAYRKVLLRKNGNKGAGGHSYGADL

There are three ways to do this,

If you click on this link Your NCBI examples. Its done.
You would then just click "download"
Click "next", click "next"

The sequence then downloads.
Essentially all I did was in the website click on the left hand navigation bar and clicked "accession number" and type in MZ039713 MZ039682 OV522182.

The second way is to go to https://ncbi.nlm.nih.gov/nuccore
Type type in MZ039713 MZ039682 OV522182
Click on the "Summary" and click on "Fasta"
Click on "send to" and click on "File"

I'm not certain is point 3 immediately above is required.
The final approach is to use efetch to automate the download, but that requires downloading the NCBI's efetch suite, which is then flipped through on a bash loop.
efetch -db nuccore -id MZ039682 -format fasta

The link for looking to automate is:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?SeqType_s=Nucleotide&VirusLineage_ss=SARS-CoV-2,%20taxid:2697049&ids=MZ039713%20MZ039682%20OV522182%20
Just to highlight
MZ039713%20MZ039682%20OV522182%20
The API is here, but I don't use them.

Thanks! I clarified that I would like the download link to be programmatic = API, not needing to do any clicking in a web app.
So do I understand it correctly that there really is no simple API to download these fastas? I need to download some efetch suite to do programmatic = non-interactive downloads?

source: Stack Exchange

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

sort: votes


text:
BWA mem | samtools view: Intermittent parsing error
Update
The issue was that bwa was running out of memory and failing, but that error wasn't floating to the top (see @Steve's answer, below). I was getting an error from samtools. I should have been able to figure that out, but hopefully this post saves someone else the trouble.
Original post
I'm experiencing a strange intermittent parsing error when piping bwa mem to samtools view. It only seems to happen when the command runs within a Nextflow pipeline (though it seems unlikely that Nextflow has anything to do with it). If I cd to the directory and run bash .command.sh, it runs to completion without issue.
Nextflow process
Here's the Nextflow process I'm running (code credit to @Steve, though modified to handle .bam or .cram). The portion that's running in this case
process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    """
    if [ "$params.output_format" = "bam" ]; then
    bwa mem \\
        -p \\
        -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
        "${params.align_to_ref}" \\
            "${fastq}" |
    samtools view \\
        -1 \\
        -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    else
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
        -M \\
        -C \\
        -H <(grep "^@RG" "${header}") \\
            "${params.align_to_ref}" \\
        "${fastq}" |
        samtools view \\
            -C \\
            -T params.align_to_ref \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
    fi
    """
}

Exact error
Here's an example of the intermittent errors I'm getting. I would have assumed that some read (e.g., read # 2_102_429 in the .fastq) is malformed, but as I mentioned, running bash .command.sh in the working directory for the failed run finishes without issue.
Command error:
  [M::mem_pestat] skip orientation FF
  [M::mem_process_seqs] Processed 700000 reads in 109.002 CPU sec, 16.411 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (24, 301945, 4, 10)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (118, 180, 399)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 961)
  [M::mem_pestat] mean and std.dev: (179.38, 124.36)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 1242)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (247, 319, 420)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 766)
  [M::mem_pestat] mean and std.dev: (342.90, 127.88)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 939)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (326, 764, 1452)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 3704)
  [M::mem_pestat] mean and std.dev: (893.30, 530.98)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 4830)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 97.823 CPU sec, 14.292 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (14, 196658, 9, 15)
  [M::mem_pestat] analyzing insert size distribution for orientation FF...
  [M::mem_pestat] (25, 50, 75) percentile: (206, 422, 1021)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2651)
  [M::mem_pestat] mean and std.dev: (506.92, 458.91)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3466)
  [M::mem_pestat] analyzing insert size distribution for orientation FR...
  [M::mem_pestat] (25, 50, 75) percentile: (248, 322, 429)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 791)
  [M::mem_pestat] mean and std.dev: (345.17, 129.22)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 972)
  [M::mem_pestat] skip orientation RF as there are not enough pairs
  [M::mem_pestat] analyzing insert size distribution for orientation RR...
  [M::mem_pestat] (25, 50, 75) percentile: (485, 777, 1315)
  [M::mem_pestat] low and high boundaries for computing mean and std.dev: (1, 2975)
  [M::mem_pestat] mean and std.dev: (961.60, 569.27)
  [M::mem_pestat] low and high boundaries for proper pairs: (1, 3805)
  [M::mem_pestat] skip orientation FF
  [M::mem_pestat] skip orientation RR
  [M::mem_process_seqs] Processed 700000 reads in 275.897 CPU sec, 39.910 real sec
  [M::process] 0 single-end sequences; 700000 paired-end sequences
  [M::process] read 700000 sequences (70000000 bp)...
  [W::sam_read1_sam] Parse error at line 2102429
  samtools view: error reading file "-"

A previous sample complained about an @SQ header line missing the LN: tag. It's like the pipe gets interrupted and lines get truncated, or something.
  [E::sam_hrecs_update_hashes] Header includes @SQ line "KI270330" with no LN: tag
  samtools view: failed to add PG line to the header

Anyone seen this before?

I haven't seen that particular error before, but I'd first make sure your pipeline errors aren't being masked by your shell. Consider adding the following to your nextflow.config:
process {

  shell = [ '/bin/bash', '-euo', 'pipefail' ]
}

IIRC, the default shell (as provided by Nextflow) does not include the pipefail option for pipelines.

The only other thing I can think of is to make sure your reference FASTA (and BWA index files) are localized in the workDir. You might find the intermittent (filesystem?) errors maybe go away even if you are staging using symlinks.
params.ref_fasta = './bwa_index/ref.fasta'

align_to_ref = Channel
    .fromPath( [ params.ref_fasta, "${params.ref_fasta}.*" ] )
    .collect()

process bwa_mem_proc {

    tag { "${sample_name}:${fastq.name}" }

    label 'bwa_mem_proc'

    input:
    tuple val(sample_name), path(fastq), path(header)
    path align_to_ref

    output:
    tuple val(sample_name), path("${fastq.baseName}.unsorted.mini.bam")

    script:
    def task_cpus = task.cpus > 1 ? task.cpus - 1 : task.cpus

    if( params.output_format == "bam" ) {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -1 \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    } else {
        """
        bwa mem \\
            -p \\
            -t ${task_cpus} \\
            -M \\
            -C \\
            -H <(grep "^@RG" "${header}") \\
            "${align_to_ref.first()}" \\
            "${fastq}" |
        samtools view \\
            -C \\
            -T "${align_to_ref.first()}" \\
            -o "${fastq.baseName}.unsorted.mini.bam" \\
            -
        """
    }
}

Note that in the first case (where params.output_format == "bam"), the reference FASTA file itself is not required by BWA MEM. Only the BWA index files are needed here.

Indeed! You were right again. Turns out bwa was running out of memory but that error wasn't floating to the top.

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


text:
Secondary structure assignment tool
Is there any CA-based online tool/web-server available for secondary structure assignment of proteins except DSSP or STRIDE?

DSSP is a database of secondary structure assignments (and much more) for all protein entries in the Protein Data Bank (PDB)
STRIDE Protein secondary structure assignment with stride, Basic assignment, Visual assignment, Contact map, Ramachandran Plot


text:
How to summarize multiple exon copy numbers into copy number of the corresponding gene
I have a matrix, sample by exon, containing a copy number value for each pair (sample, exon). I would like to generate a second matrix, sample by gene, where the copy number of the exons is consolidated in a single value per gene.
Due to the nature of copy number alteration, I am not sure if this is meaningful, or feasible (e.g., how can it be summarized, in the case a gene is interested by 2 alterations with different copy numbers?).
Is there a common approach that is used in this case?
Edit
A possible example is the following:
         gene1_exon1  gene1_exon2  gene1_exon3 ... geneN_exonM
sample1            2            1            1               0
sample2            2            2            1               0
...
sampleJ            3            2            2               2

A possible summarized matrix (with majority voting - not sure if this would be the best approach):
         gene1 ... geneN
sample1      1         0
sample2      2         0
...
sampleJ      2         2

I believe this is more an issue related to the biological interpretation than the actual implementation (i.e., I am already implementing a function in pandas). My concern is on the interpretation, since my prototype select as gene copy number the majority copy number from its exons. However, this has some caveats (e.g., what if you have a GAIN on an exon, and a LOSS on a second exon, for the same gene - will you call GAIN or LOSS?). The big picture is that I want to compare two samples CNAs at the gene level, but I could not find this task discussed online.
The main objective of this summarization is to compute gene level concordance between samples. Exon-level or base-level concordance is computed in parallel, in order to have different CNA granularities.
However, one might not use the majority voting proposed and, e.g., call a CNA on a gene if the alteration is found in at least p% of exons (e.g., 30% of exons) - this is just an example, I am not suggesting to use this approach.

You might be interested in the approach we took when trying to find CNVs that overlap known, reported ones. have a look at our documentation here: https://varsome.com/about/resources/sv-implementation/. I'll try and get the relevant bits into an answer if I get the chance.

Thank you! I will look into it

Take it with a pinch of salt, of course: that was not peer reviewed and was just the best system we could come up with based on our own best judgment. That said, I do think it is a reasonable approach, I am just aware that I have been wrong before! :)


text:
How to subset a GRanges object based on a specific genomic window of interest?
Goal: To keep only rows of a GRanges object if their start coordinate is within 1 Mb of the start coordinate of a gene of interest.
I am aware of some functions that come close:

GenomicRanges::nearestKNeighbors() (and related distanceToNearest(), nearest(),
precede(), follow())
GenomicRanges::subsetByOverlaps()
GenomicRanges::window()

And I am aware of the UNIX command-line 'Bedtools' tools: window, closest, and intersect, however I wish to accomplish this task using R only and the above functions have not been doing exactly as I would like.
Here is an example GRanges dataset:
> gr <- GRanges(
        seqnames = Rle(c("chr1", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
        ranges = IRanges(start = seq(300000, 2000000, length=10), end = seq(700000, 6000000, length=10), names = head(letters, 10)),
        strand = Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
        score = 1:10, GC = seq(1, 0, length=10), GeneSymbol = c(paste("gene", seq(1:10), sep="")))

> gr
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  a     chr1   300000-700000      - |         1  1.000000       gene1
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  e     chr1 1055555-3055555      * |         5  0.555556       gene5
  f     chr1 1244444-3644444      + |         6  0.444444       gene6
  g     chr3 1433333-4233333      + |         7  0.333333       gene7
  h     chr3 1622222-4822222      + |         8  0.222222       gene8
  i     chr3 1811111-5411111      - |         9  0.111111       gene9
  j     chr3 2000000-6000000      - |        10  0.000000      gene10
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Ideally, by specifying the name (e.g., "gene3") or start coordinate (e.g., chr2 AND 677777) of the gene of interest AND a window size (e.g., 1000000 or 1 Mb upstream AND downstream) to a function, I hope to get an output like below:
GRanges object with 3 ranges and 3 metadata columns:
    seqnames          ranges strand |     score        GC  GeneSymbol
       <Rle>       <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2  488888-1288888      + |         2  0.888889       gene2
  c     chr2  677777-1877777      + |         3  0.777778       gene3
  d     chr2  866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Notice that only three rows from the original GRanges object remain in this example output because the start coordinates for all other genes were either (i) on different chromosomes or (ii) were more than 1 Mb away from the start coordinate of gene3.

I have tried various methods to isolate the rows within 1 Mb of my gene of interest with no luck:
> gr1 <- gr[grep("gene3", gr$GeneSymbol),]
> grOut <- subsetByOverlaps(gr, gr1, maxgap = 1000000, invert=T)

> grOut <- gr[mcols(gr)$GeneSymbol == "gene3" & start(gr) < 677777-1000000 & end(gr) > 1877777+1000000]

> grOut <- window(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)))
Error: 'start' must be a single number or NA

> grOut <- gr[IRanges(start=677777, end=1877777)]
Error: subscript contains out-of-bounds ranges

> grOut <- nearest(gr, GRanges("chr2", IRanges(start=677777-1000000, end=1877777+1000000)), select="all", ignore.strand=T)

The last attempt recorded above was likely the closest I've gotten to my desired output by storing the index of rows in the GRanges object falling within the range specified, but it presents this index in the form of a SortedByQueryHits object, not a GRanges object as I require. It seems simple enough a task, but I haven't solved it yet; any thoughts or solutions?

For those wishing to find/verify the coordinates for a gene of interest (GOI) in their GRanges object when only the gene symbol/ID/name is known: use gr[grep("<GOI>", names(gr)),] when column names have been set as gene names ... use gr[grep("<GOI>", rownames(gr)),] when row names have been set as gene names ... use gr[grep("<GOI>", mcols(gr)$<geneNameColumn>),] when a metadata column contains gene names.

Specifying the range of interest (+/- 1 Mb of start coordinate) as a new GRanges object and then subsetting the original GRanges object based on the overlap of genes with the query works:
library("GenomicRanges")

> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]
> range <- 1000000

> q <- GRanges(seqnames=chrPos,
               ranges=IRanges(start=startPos-range, end=startPos+range))
> grOut <- subsetByOverlaps(gr, q)
> grOut
GRanges object with 3 ranges and 3 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character>
  b     chr2 488888-1288888      + |         2  0.888889       gene2
  c     chr2 677777-1877777      + |         3  0.777778       gene3
  d     chr2 866666-2466666      * |         4  0.666667       gene4
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

Credit to @zx8754 for inspiring this solution.

If preferred, this can be achieved with base R functions:
# First store the chromosome and start position of your gene of interest in variables
> chrPos <- as.character(seqnames(gr))[gr$GeneSymbol=="gene3"]
> startPos <- start(ranges(gr))[gr$GeneSymbol=="gene3"]

# Then add a column to the GRanges object with distances from the start position of each
# gene in the GRanges object to the start position of your gene of interest
> gr$distance <- start(ranges(gr))-startPos

# Finally subset the GRanges object based on your desired window of interest (must be
# symmetrical —— e.g., +/- 1 Mb —— for the following to work)
> grOut <- gr[seqnames(gr) %in% chrPos & (abs(gr$distance) <= 1000000)]
> grOut
GRanges object with 3 ranges and 4 metadata columns:
    seqnames         ranges strand |     score        GC  GeneSymbol  distance
       <Rle>      <IRanges>  <Rle> | <integer> <numeric> <character> <integer>
  b     chr2 488888-1288888      + |         2  0.888889       gene2   -188889
  c     chr2 677777-1877777      + |         3  0.777778       gene3         0
  d     chr2 866666-2466666      * |         4  0.666667       gene4    188889
  -------
  seqinfo: 3 sequences from an unspecified genome; no seqlengths

This strategy provides the added benefit of adding a 'distance-from-gene' (gr$distance) column to the dataset which can be useful for visualizations.


text:
What is the best way to compare treated/control samples?
I have two populations of mice, say RED and GREEN; there are 5 pairs treatment/control in each.
The table is the following:
gene_name  p-value(treatVSctrl) RED1_ctrl RED1_treat RED2_ctrl RED2_treat ... GRN1_ctrl GRN1_treat etc...
gene1      1.3                  4         10         5         10             6.2       9.75
gene2      2.8                  3         6          6         12             3         6
gene3      4.4                  5         12         4         8              7         10

I need to compare two groups to have an idea of the expression of different genes after treatment.
I also need to find a difference between RED and GREEN group in terms of genes expression
The thing is I cannot compare mice one by one just measuring the fold change, I need to consider p-value. What is the most correct way to perform this analysis?

Its an easy comparison but I was confused by the header p-values. I can't help there needs to be raw quantification. I am not aware of a stat that can take an unknown 'p-value' and compare to other p-values at a given level of significance

Thank you. Do you mean there must be a p-value for each pair?

No the header in the first column says " p-value(treatVSctrl..." . I don't understand the quantification, but the column header is stating this is a p-value.

Agree the "consider the p-value" instruction is confusing. The first column clearly doesn't represent p-values, as p-value are probabilities and must be in the range [0, 1].

Is this RNA-seq? How many genes are there in total and how many samples per group? The suggested Wilcox is not very powerful so it would only work if your group size is sufficient, and if this is RNA-seq then you sould use established software such as DESeq2 or edgeR which takes care of the stats and all the necessary internals. Are these data normalized already, what is the experiment?

@NuclearHoagie it could represent a percentage, i.e. (0 to 1.0) * 100 , e.g. this could be a pairwise T-test. Alot of values are hovering around the 5% mark. Where there are mixed results it would be a reason for the OP to ask how to move forward.

I don't really understand the table, particularly the 'p-value' in the header. More generally I don't really understand the experiment nor the data.
As a general rule the Wilcoxon would be a reasonable guess of the test using raw data. This is non-parametric version of the T-test without the rigid assumptions of parametrics. Control mice go in one column and subjects in the other. Scipy.stats has the test, its certainly in R (which I don't really use).
print (scipy.stats.wilcoxon(x = [1,2,3,4,5], y=[1,3,5,6,7], zero_method='wilcox', correction=False, alternative='two-sided', mode='auto', *, axis=0, nan_policy='propagate'))

There you go, filled out with some example data.


text:
GATK variant filtration using "SelectVariants" and use of JEXL queries
How to do variants selection in some corner cases using GATK and JEXL expressions?
I am following the guidelines given in this links for variant selection for some specific cases are not as much detailed.

https://gatk.broadinstitute.org/hc/en-us/articles/360036430432-SelectVariants
https://gatk.broadinstitute.org/hc/en-us/articles/360035891011-JEXL-filtering-expressions

I just had to do some filtering of my own data and tested several methods (to see what the outcome would be). I am sharing those here for others and myself, if it is ever needed again.

I am following the guidelines given in this links for variant selection for some specific cases that are not as much detailed in the GATK site.

https://gatk.broadinstitute.org/hc/en-us/articles/360036430432-SelectVariants
https://gatk.broadinstitute.org/hc/en-us/articles/360035891011-JEXL-filtering-expressions

I just had to do some filtering of my own data and tested several methods (to see what the outcome would be). I am sharing those here for others and myself, if it is ever needed again.
Tutorials .... !!!
- method to call the sites that are all homozygous Variants (1/1, 2/2 genotype)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() == 6'
# outcome: no sites will have ambigous (./.) and/or homRef (0/0). All other homVar are allowed

- in the above script if we want the site where at least 4 samples are homozygous ref or variants
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() >= 4'
# outcome: atleast 4 samples are homVar, rest can be any genotype (./., 0/0, 0/1, etc.)

- at least 1 sample is homVar (1/1, 2/2 or 3/3)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() != 0'

- all the samples are homRef (0/0 genotypes)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomRefCount() == 6'

Note: 6 chr = 3 samples (when samples are diploid)
- sites that are all HomVar or nocalls (./.)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() == vc.getCalledChrCount()/2'
# outcome: no 0/0, 0/1
# Note: samples with GT = ./. are not counted as called chr

Useful after samples are sub-sampled from multisample vcfs
- sites where no genotypes were called for all the samples in the selected vcf
-  or say, select vcf sites that are not called in any sample (i.e don't have any genotypes, only contain ./.)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() == 0'
# outcome = vcf sites (lines) that are all no calls (GT = ./.) on all samples are selected

- samples are either homRef (0/0) or no call (./.)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomRefCount() == vc.getCalledChrCount()/2'
# outcome: vcf sites where all the samples are GT = ./. (no call) or GT = 0/0

- site where 6 chromosomes (exact) or 3 samples were called for diploid genome
java -jar /home/everestial007/GenomeAnalysisTK-3.7/GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() == 6'
# Note of Caution: 6 chr = 3 samples (when samples are diploid)

- site where at least 3 samples (diploid) are called
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() >=6'

- if I have 6-sample vcf
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getCalledChrCount() ==12'
# 12 haploid chr = 6 samples (diploid)
# outcome: sites that were called in all samples, there will be no GT = ./.

- select vcf sites where no sample has hetVar (no GT = 0/1, 0/2), only GT = ./., 0/0, 1/1, 2/2
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHetCount() == 0'

- select vcf sites where no sample has hetVar (no Gt = 0/1, 0/2), and at least one sample is HomVar (GT = 1/1)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHetCount() == 0 && vc.getHomVarCount() >= 1'

- site where at least 1 sample is homVar (GT = 1/1 and/or 2/2 etc.)
java -jar GenomeAnalysisTK.jar -T SelectVariants -R lyrata_genome.fa -V MY.phased_variants.Final.vcf -o My.filtered_01.vcf -select 'vc.getHomVarCount() >= 1'

- where a site is only homRef for a particular sample and is only SNP
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getGenotype("ms01e").isHomRef()' -O test_variants.ms01e.Chr_2_3.vcf

- where a site is only homVar for a particular sample and is only SNP
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getGenotype("ms01e").isHomVar()' -O test_variants.ms01e.Chr_2_3.vcf

- where a site is not homRef for a particular sample and is only SNP. This will also select empty sites
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select '! vc.getGenotype("ms01e").isHomRef()' -O test_variants.ms01e.Chr_2_3.vcf

NOTE: following genotype call methods won't work:
# to get only heterozygous sites
-select 'vc.getGenotype("ms01e").isHetVar()

# to get any variant sites
-select 'vc.getGenotype("ms01e").isVar()

# won't work if getGenotype method is combined
-select 'vc.getGenotype("ms01e").isHomVar() || vc.getGenotype("ms01e".isHomRef()'

So, I figured out the following methods where count helps with selecting particular site
- where a site is a heterozygous variant site for a particular sample
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select 'vc.getHetCount() == 1'  -O test_variants.ms01e.vcf

NOTE: 'vc.getGenotype("ms01e").isHetVar()' wont work so generated the above expression.

- where a site is a variant site (either homVar or hetVar) for a particular sample and is only SNP
java -jar GenomeAnalysisTK.jar SelectVariants -R lyrata_genome.fa -V phased.MySpF1.vcf -sn ms01e -select-type SNP -select 'vc.getHetCount() == 1 || vc.getHomVarCount() == 1'  -O test_variants.ms01e.vcf

NOTE: 'vc.getGenotype("NA12878").isHomRef()' || vc.getGenotype("NA12878").isHomVar()' wont work together, so generated this above expression.

Other updates to follow.
Thanks,

Restored because this is an answer to the question [i.e. a user has answered their own question].


text:
tRNAdb insertion symbols. To what they correspond?
I have to manually curate some tRNA sequences in order to have their original sequence without modifications. I am getting these sequences from tRNAdb. In these sequences, I found often "_" symbols.
They provide this list of symbols and "_" corresponds to insertion. My question is the insertion of what? Are they referring to a reference sequence? In the list they say that in the footnotes they provide more info. However, there are no footnotes.
Have you ever had experience with tRNAdb?

I have BLASTed one of the tRNA sequences. Lactobacillus acidophilus, Ile, that has insertion sites.

GGGCCTATAGCTCAGCTGGTTTAGAGC_CACGCCTGATAAGCGTG_GTCGATGGTTCAAGTCCATTTAGGCCCA

It looks like insertion refers to sites that exist in reference DNA but missing from tRNA sequence. But I am not a tRNA expert, it's my interpretation. See the image.

I post the workaround that I eventually found for this problem as it could serve as a reference for the exploration of tRNAdb bioinformatic resource.

If you look at the picture, tRNAdb provides comments about the insertion symbols. Therefore, for the specified position, nucleotides are defined. The original tRNA sequences can be reconstructed from this information.
Hope that this problem could be of help to the community!


text:
How can I list all amino acids in a protein using PyMOL?
What is the command in PyMOL for listing all the amino acids in a specific protein, say, 1a62.pdb?

At the command, change the current directory to your target directory:
cd my-directory-path-containing-my-files

Load the file:
load 1a62.pdb

Then you can either print to screen or save to a file:
print(cmd.get_fastastr('all'))

save 1a62.fasta

A small caveat for get_fastastr is that the fasta sequence returned does not include gaps in the numbering, therefore, a multiline command using [int(atom.resi) for atom in cmd.get_model('name CA').atom] may be required in some cases.


text:
Nextflow scheduler job doesn't end upon successful workflow completion
Bug fixed
This issue appears to have been fixed in latest Nextflow: https://github.com/nextflow-io/nextflow/issues/2693

I've been using Nextflow for various workflows on a cluster using the SLURM scheduler. I've noticed that the 'parent' job does not shut down when the Nextflow workflow finishes successfully. I can cancel the job, but then Nextflow doesn't get a chance to clean things up and also doesn't write out the trace, report, etc. The parent job does end properly if the workflow ends with an error.
Question: Why doesn't the parent Nextflow job end upon successful completion?
I saw a post on the Nextflow GitHub where the original poster describes various acrobatics to get the parent job to exit cleanly, but this is something that Nextflow should handle.
Here's my SLURM submission script that I submit using sbatch run_workflow.slurm
run_workflow.slurm
#!/bin/bash
#SBATCH --time=60:15:00                             # Time limit for the job (REQUIRED).
#SBATCH --job-name=NF_Parent_Identify_D&C_regions   # Job name
#SBATCH --ntasks=1                                  # Number of cores for the job. Same as SBATCH
#SBATCH --mem=5G                                    # Total memory requested
#SBATCH --partition=<partition>                     # Partition/queue to run the job in. (REQUIRED
#SBATCH -e slurm-%j.err                             # Error file for this job.
#SBATCH -o slurm-%j.out                             # Output file for this job.
#SBATCH -A <account>                                # Project allocation account name (REQUIRED)

module load ccs/java/jdk1.8.0_202
nextflow run identify_dark_and_camouflaged_regions.nf \
    -with-report \
    -with-trace \
    -with-dag idcr-workflow.png \
    -with-timeline \
    -resume

Update 1
I'm adding example Nextflow output from a different workflow (I don't have the SLURM output from the one I mention above) where the workflow has completed, but the parent job is still running. In this case, there were errors, but I added errorStrategy 'ignore' to all of the processes listed:
executor >  slurm (876)
[06/c5672b] process > samtools_csort_proc (A-CUHS... [100%] 304 of 304, faile...
[45/26a046] process > Lsamtools_csort_proc (A-CUH... [100%] 304 of 304, faile...
[b5/0f2d7c] process > sambamba_csort_proc (A-CUHS... [100%] 304 of 304, faile...

So, regardless of whether the workflow completed successfully, Nextflow has nothing left to do, but it's still running.
Update 2 (adding partial .nextflow.log file)
Per @Steve's point, I checked the .nextflow.log file and am including part of the file, here. In this case, I don't think Nextflow has any tasks related to publishDir because I'm not using that in this workflow. The job that's still running is from the Nextflow workflow in my GitHub repository for my answer posted in my related question Samtools sort: most efficient memory and thread settings for many samples on a cluster. Whether it's reproducible, I haven't checked. Not sure how common this issue is.
In either case, my interpretation of the .nextflow.log file below is that Nextflow thinks it's hung, too. I just don't know why, or how to exit cleanly.
Partial view of .nextflow.log
...

Feb-28 05:29:04.607 [Task monitor] DEBUG n.processor.TaskPollingMonitor - !! executor slurm > tasks to be completed: 1 -- submitted tasks are shown below
~> TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: RUNNING; exit: -; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: -; ]
Feb-28 05:30:54.613 [Task monitor] DEBUG n.processor.TaskPollingMonitor - Task completed > TaskHandler[jobId: 662274; id: 876; name: sambamba_csort_proc (A-CUHS-CU003128-BL-COL-49696BL1.unsorted;CPUs:2;MEM_PER_THREAD:100GB); status: COMPLETED; exit: 0; error: -; workDir: /path/to/Samtools_sort_optimization_test/SORT_COMPARISON/work/b5/0f2d7c5d7f81a6e62dedb08bf6694a started: 1646038174334; exited: 2022-02-28T10:30:53.516015Z; ]
Feb-28 05:34:04.641 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:39:04.647 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:44:04.653 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:49:04.658 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled
Feb-28 05:54:04.662 [Task monitor] DEBUG n.processor.TaskPollingMonitor - No more task to compute -- Execution may be stalled

< above line repeats hundreds of times >

...

Can you please create an answer with your accepted solution, so that the StackExchange algorithms don't confuse this with an unanswered question?

How do you know your workflow actually finishes successfully? One way to check the workflow termination status is to use the workflow.onComplete and workflow.onError event handlers. Consider resuming your workflow with something like the following appended to your script:
workflow.onComplete {

  log.info "Workflow completed at: ${workflow.complete}"
  log.info "Time taken: ${workflow.duration}"
  log.info "Execution status: ${workflow.success ? 'success' : 'failed'}"
}

workflow.onError {

  log.info "Execution halted: ${workflow.errorMessage}"
}

Note that the workflow event handlers may also be defined inside your nextflow.config if preferred.

Thanks, @Steve. I added an example to my post. Seems that one way or another, the example I show above should have completed because Nextflow has nothing left to do.

@MarkEbbert It's hard to say what the problem is here exactly without the nextflow.log or a reproducible example. But just because all of the tasks have completed (regardless of the exit status) doesn't necessarily mean Nextflow has nothing left to do. For example, the publishDir directive works asynchronously and can be used (as you know) to 'copy' task outputs. If these are large and there are many of these, Nextflow will hang around copying files from the queue until well after the last task has completed.

Thanks, @Steve. I just added part of the .nextflow.log file. Looks like Nextflow thinks it's hung, too? I haven't tested whether this is reproducible, by this is the code for my related question. I included a link to it above. This one may just remain a mystery.

I can't edit my above comment, but should read "..., but this is from the code..."

@MarkEbbert Correct - this looks like a bug. It might have something to do with the errorStrategy set to 'ignore' for all processes, potentially. Basically, you need to hit this condition to break out of the poll loop. From the partial log you've posted, it looks like the runningQueue is empty, but I don't know about the pendingQueue.size() or if session.isTerminated().

@MarkEbbert You might like to open an issue with the full log attached.

Thanks, @Steve. Submitted: https://github.com/nextflow-io/nextflow/issues/2693. I'll report back if I learn something.


text:
Replacing column 1 in file1 with matching ID's in file 2
I have downloaded two files:
File1.csv from:
https://depmap.org/portal/download/data_slicer/download?file_path=20220227%2F3c335b7d-c18e-4788-b2e6-512d71f0a339%2Fexport.csv&name=Metabolomics.csv
File 1 looks like
,2-aminoadipate,3-phosphoglycerate,alpha-glycerophosphate
ACH-000698,6.1127272,6.034197900000001,5.8968963
ACH-000489,5.5774126,5.7270453,5.1114679,6.07325

File2.csv from:
enter link description here
https://cog.sanger.ac.uk/cmp/download/model_list_20220205.csv
File2 looks like:
 model_id,model_name,synonyms,model_type
 SIDM01774,PK-59,,Cell Line,Adherent
 SIDM00192,SNU-1033,,Cell Line,Adherent

I want to replace ACH-* ids with SIDM* ids in File1.csv. The ACH-* id's for SIDM* ids are given in the "BROAD_ID" (col.no =40) column in the File2.csv
Desired output
,2-aminoadipate,3-phosphoglycerate,alpha-glycerophosphate
SIDM00523,6.1127272,6.034197900000001,5.8968963
SIDM00835,5.5774126,5.7270453,5.1114679,6.07325

I want to do this with python script. I tried some initial steps. Please forgive any grave mistake, as I am a beginner in Python
 import pandas as pd

 df = pd.read_csv("File1.csv", index_col=0).T

 print(df)

 models = pd.read_csv("File2.csv", index_col="model_id", usecols=["model_id", "BROAD_ID"]).rename(columns={"model_id": "Model"})

 df = df.join(models).rename_axis(index="Model Name").reset_index()

 print(df)

 df2 = df.set_index("Model").drop(columns="Model Name")

 df2.to_csv("Meta.tsv", sep="\t")

 df2.index = df2.index + ".Demeter"
 df2.T.rename_axis(index="Gene").to_csv("Meta_2.tsv", sep="\t")

I was not sure about the exact output (columns) you desired. The following code should inner merge both data frames. It writes a copy of 'model_list' column as a first column. You could remove any unwanted columns. I am sure there are more elegant ways of doing it. Hope it helps. It's totally different approach than yours.
import pandas as pd

metabolics = pd.read_csv("Metabolomics.csv", sep=",", header=0)
model_list = pd.read_csv("model_list_20220205.csv", sep=",", header=0)
metabolics=metabolics.rename(columns={'Unnamed: 0':'ACH_number'}) #rename the unnamed column.
merged_df = pd.merge(metabolics, model_list, left_on='ACH_number', right_on='BROAD_ID', how='inner') #join two files

#insert model_id column at the beginning
#You will have dupicae 'model_id' column
merged_df.insert(0, "model_id", merged_df['model_id'], allow_duplicates=True)
merged_df.to_csv('results.csv', sep=",", index=False)
print(merged_df.shape)

This looks about right

I agree with @M__ 's comment below. You would need to test the code make sure that merge on keys/columns as indicated give the desired output.

My pandas is rusty. The basic issue is that joint in SQL is merge in pandas. This might be why the OP used joint (except its pandas).  I was always told SQL joint merges via a common key and I carried that through into pandas.

I am no expert either. But I will look into this as I start using it recently more often.

Hi. I used the above script. There were two issues:

I am getting a result file with 973 rows, while the original file have 929 rows.
I want only rows where a SIDM* code is available for ACH-. Also in my desired output I don' t want the ACH- code column.

Code should be do the right thing. If you inner merge, your results file will have 923 lines excluding headers. Now, you file has 928 record. You are missing the following 5 records from results.> ACH-000338, ACH-000467, ACH-000710, ACH-000833, ACH-001063. The reason for is that one moodel_list file has irregular entries. For example if you look at the col-40, you will see that actual cell entry is "ACH-000338;ACH-000338", not "ACH-000338" thus pandas cannot merge these.

Either you live losing 5 records, or you need write another script that considers this kind of anomalies or update the tables. Hope this helps.

@Supertech: Thank you so much for getting back. The script is correctly getting the 923 matching id's. But when I saw the first column with awk. There were extra entries in the first column from the random columns of model_list file

The basic strategy of merge 'left' (which I think is default) is the correct thing to do. However, I would be concerned by lack of a key between the two data sets and basically the merging will take place on the row order. Within merge is on= then the shared key, in the code by @supertech there is left-on='ACH_number', I dunno because ACH_number does not appear to be shared. The code needs testing.
If you assign an index it will attempt to merge on the index (which I think will therefore fail). What I think you want is the row order 1,2,3,4 in each dataframe to be the key that is merged. Therefore not specifying key might be the right thing to do, what can happen is it will then use the row order as the left-key. I cannot remember.
The insert is correct to reorder the columns to get the desired output.
Finally you want,
del merged_df['column_name']

for each column you don't want. There is a way to loop through this, but the above keeps things simple. You then do the standard csv dump.
Note, see comments on drop by @GOATNine, which strongly preferred.
Setting and unsetting indexes is key to manipulating pandas but it is within a specific contex.
Summary the command merge and insert are correct with the addition of del however, the code needs testing regarding merge left-on and assigning the index. It is much easier to print the output and visualise it.

Thank you so much for explaining it so well

I do not know the reason why wget fails. Maybe someone else could help. I have clicked on the link on my Pc and it was sufficient to download. wget also failed for me.

If you're using pandas to drop the columns, you can just use df = df.drop(columns=['column1','column2'....]), or drop the assignment and use the option 'inplace=true' instead. This allows you to drop as many columns as you like at once, by making a string list of the column names.

Yes, I remember now. This is much better

Could anyone please update the code. I am a begginer in python and I am not getting where should I include all the steps suggested. I want file1 with all the columns except column 1 (ACH-* ids) replaced by SIDM* id's.

merge, insert, drop. I might test it if there are example tables (not links), or I was currently working with pandas. The info you supplied isn't clear enough for me so I can easily transfer it an IDE

@Megha see my comment under my answer. Your model_ file has anomalies in some of the cells therefore inner merge fails for 5 of the records.


text:
SNPs in LD in which populations?
I have 3 SNPs that are in LD in 'several' populations here:
https://www.ebi.ac.uk/gwas/variants/rs13078854#ld_panel_label
It's a pain to look at all (30) populations one by one to manually check if the 3 SNPs are in LD in a given population or not.
Is there a way to ask this question more systematically? e.g. in which common imputation populations are a given set of SNPs in LD (above a certain threshold)?

There is a risk of a type 2 error.
My advice would be to use Structure if it is still available given the situation. Its not great and it is old and its a GUI and I don't know whether it has been updated (I'm not exactly selling this am I?), but it will do the job of navigating through a type 2 error.
What you do is provide the number of populations a priori and it will assign a Bayesian estimate of HW equilibrium for that number of populations. You manually input 1,2,3,4,.... 30 run and wait. The point at which the joint marginal plateaus is the minimum number of populations which do not exhibit LD (thus if two of your populations exhibit LD it collapse them into one). Thus if your answer is '20' it give some confidence it is a real result. If it is like 10 (optimistically) the conclusion would be in pretty reasonable shape.
There are other ways of looking at this, but that would be my starting point.
The classic calculation is LDHat, here https://github.com/auton1/LDhat, but the last update was 2015 and in this particular circumstance I'm not sure.

I don't understand.

Sorry if this has not been helpful. I would therefore suggest a data pipeline around what appears the current standard which is plink https://www.cog-genomics.org/plink/1.9/ld . To understand LD stats a background in population genetics and stats is really helpful. The answer above is as good as I could do without writing a data pipeline.

We can chat about it, but the solution I proposed was trying to over the type 2 error in the data and I'd need to explain that in depth because its fundamental to how LD works. I'm not up to date with plink, I left LD prior to plink, but I do understand how the numbers work in LD.

How would you like to chat?

I forget how you do invites and stuff, but I think you have to visit the bioinformatics chat and that activates an invite.


text:
samtools view command not found error
When I tried to use samtools to split a bam file based on different chromosomes, I used this command:
samtools view input.bam -b chr21 | chr21.bam

However, I get error messages like this:
-bash: chr21.bam: command not found
[W::hts_idx_load3] The index file is older than the data file: input.bam.bai

How could I solve this problem? Thanks in advance~~~
The version I am using is:
Program: samtools (Tools for alignments in the SAM format)
Version: 1.10 (using htslib 1.10.2)

Try re-indexing your bam file:
samtools index input.bam

Or update the timestamp via touch index.bam.bai but I think re-indexing is safer. You might want to save any original index file first, e.g. via mv input.bam.bai old_input.bam.bai or similar
Also I think you meant > chr21.bam instead of piping via | chr21.bam if you want to save the output file

Cool! I think the problem is because I use | chr21.bam.


text:
bcftools mpileup "Failed to open the file" questions
When I trying to use bcftools mpileup to convert bam to vcf files,
I try to use:
bcftools mpileup -Ob -o resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai input.bam

However, I got this error message:
[E::fai_build3_core] Failed to open the file resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai

Or, I would get another error message when I use another ref file:
[E::fai_build_core] Format error, unexpected "c" at line 1

Does anyone know how to solve these problems? Thanks in advance
Also, could anyone recommend some human gene Reference? The one I am using now is from GATK pipeline (but I don't know if that is good enough).
The version I am using is:
Program: bcftools (Tools for variant calling and manipulating VCFs and BCFs)
License: GNU GPLv3+, due to use of the GNU Scientific Library
Version: 1.10.2 (using htslib 1.10.2)

What are you trying to do exactly? The -o option is to set an output file, and you seem to be trying to overwrite an existing fasta index (.fai) file. Did you mean to use -f to give the genome (which would be resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta and not resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai)?

Regarding a human genome reference, have a look at CHM13/HG002 from the Telomere-to-Telomere consortium:
https://github.com/marbl/CHM13#assembly-releases
The benefit of the T2T assembly is that it is fully-assembled across all repetitive regions, which improves the mapping accuracy. The pangenome consortium is working on additional T2T assemblies from other people, trying to capture as much human [structural] variation as possible.
https://humanpangenome.org/
The v2.0 assembly, adding chrY from HG002, has only just been released. From the looks of it, it is not yet annotated, but there's gene annotation for all other chromosomes via the v1.1 release:
https://github.com/marbl/CHM13#downloads


text:
Identify side chain atoms in BioPandas dataframe
I am currently implementing some metrics I could use for comparing two conformations of the same protein in Python.
For example, I know I could use the RMSD of all protein ATOMS in BioPandas using this:
biopandas.pdb.PandasPdb.rmsd(conformation1.df['ATOM'], conformation2.df['ATOM'])

Are there some out-of-the-box tools in Python for measuring the RMSD between the side chains of the two conformations? When trying to do this manually, the main issue is how to identify side chain atoms, which is not clear for me. I know how to get all atoms from a residue using BioPandas:
aminoacids = list(PandasPdb().read_pdb(my_pdb_filename).df['ATOM'].groupby(['chain_id', 'residue_number', 'residue_name']))

Now aminoacids[i] is a dataframe containing all the atoms of the i-th aminoacid, like this:

From here, it's not clear how to separate back bone atoms from side chain atoms. It's clear that the one with atom_number 'CA' will be the alpha-Carbon, but I am not sure how to identify the others. But once the separation is done, the side chain RMSD could be measure for example using Bio.SVDSuperimposer.SVDSuperimposer which works directly on the array of 3D coordinates of the atoms.
Thanks!

As you suggest one way of solving your problem would be by selecting all atoms that don't have backbone atoms names. In a pdb file I believe backbone atoms would be named 'CA', 'HA', 'N', 'HN' or 'H', 'C' and 'O'. Beware of the N-terminal (where the hydrogens would be named either 'H1', 'H2' & 'H3'; or 'HN1', 'HN2' and 'HN3') and the C-terminal (the other oxygen would be named 'OXT').
One way of doing it is by using the handy selections from pandas dataframes. It's probably a more elegant way than making a list of dataframes.
from biopandas.pdb import PandasPdb

bpdf_1 = PandasPdb().fetch_pdb('1t48').df['ATOM']
bpdf_2 = PandasPdb().fetch_pdb('1t49').df['ATOM']

# Function returns a list of boolean. True when it's not a backbone atom.
sel = lambda df: df['atom_name'].str.contains('[^(CA|HA|N|C|O|HN|H)]$')

# The number of atom is not the same, so I only align the first 200
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[sel][:200], bpdf_2[sel][:200])

The regex says something like "(not (CA or HA or …)) end of line". The end of line is important, otherwise any atom containing N, C, O or H in their name would be excluded (so most of them).
Note that with this method it's still possible to use the PandasPdb selectors. For example, to only consider heavy atoms (no hydrogen):
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[sel][:200],
                                bpdf_2[sel][:200],
                                s='heavy')

EDIT:
Another way to do it would be:
sidechain_rmsd = PandasPdb.rmsd(bpdf_1[:200],
                                bpdf_2[:200],
                                s='main chain',
                                invert=True)

You don't need any selection in this case, BUT this is a bad idea if you have hydrogens, as it will also align them (including the backbone ones).

If you want to align each residue one by one instead: I'm sure this gives you the tools to do it, though it would be a pretty bad idea. In this very case I'd advise to try comparing side-chains using dihedral angles (hence rotamers). You can have a look at pyrosetta for this.

Thanks for your detailed response! To check that I understood right: the H[0-9]* and HN[0-9]* are the names of the hydrogen atoms in the amino group? So I need to filter them out.

Yes, although you can get away by just filtering out H[1-3] and HN[1-3]. More hydrogen atoms would be disturbing.

My bad: H[1-3]? and HN[1-3]?, otherwise you don't cover cases where there is only H or only HN. Remember to keep the $ at the end, otherwise you'd end up matching all H.

Thanks for your help!


text:
Selecting a part of heatmap
I have a heatmap
I see a block of genes for a group of samples
How I can select a part of heatmap?
I mean is there any way to export a list of genes from rows of a heatmap in R?
I have tried which returns nothing
> library(ComplexHeatmap)

> ht = Heatmap(heat1)
> ht[1:35,1:5]

Why would you want to do that? Why not simply subset the data that generates the heatmap to the desired rows? What have you tried to do or searched to do this?

From a matrix of scaled values, how I can say which miRNAs I want before being clustered and visualized so that in the heatmap I can see a blocked of correlated miRNAs?

You can check/create clusters outside the heatmap via hclust. Keep in mind that the heatmap will change if you provided less genes/miRNAs or different samples. But probably can extract the data from ht object.

You should really read the manual before asking questions here.

You're using ComplexHeatmap, which is pretty easy to subset. Your subset operation creates a subset object of the heatmap, not of the data in the heatmap.
In case you don't want to run hclust outside of the CH, you can use row_order and column_order to re-order your data in the CH display order.
library(ComplexHeatmap)

ht = Heatmap(heat1)
r_o <- row_order(ht)
c_o <- column_order(ht)
heat1_reordered <- heat1[r_o, c_o]
heat1_reordered[1:35,1:5]

CH does say that clustering is oly "fixed" on drawing, so you may want to run this before picking the row and column orders:
pdf(NULL)
ht <- draw(ht)
dev.off()


text:
Assess ploidy estimation through SNV and CNV results
I want to assess ploidy estimation (i.e., diploid vs tetraploid) in a set of samples for which I have both results from the variant caller (i.e., SNVs) and from the copy number caller.
What are the possible approaches that can be used in this case? I know that in some in cases it might be possible to use VAF to assess ploidy estimation, but I haven't found exactly how this (possibly simple) approach is defined.
Note: I already have a caller that I am testing. My question is more related to exploratory data analysis or manual inspection of the caller's results.
VAF distribution examples might be similar to what is done in CNAqc (image below extracted from their online tutorial):

What is the underlying data quality? Would you be able to go in and visually confirm alleles in a single sample in IGV or whatever? Approaches such as here use allele frequency at heterozygous SNPs in single samples to try to address this question, so I think that's the relevant question. I have to say that I don't understand the figure or the tutorial at first glance, it appears to be looking for de novo mutations (?) and seeing their frequency.

Thank you. I have updated the question regarding the tutorial - I'm only interested in plotting the VAF distribution for the different genotypes, that is the first example I could find online to explain one of the tools for the analysis (sorry it that confused you). I will read the paper you linked, thank you.


text:
Range of values for VAF in heterozygous genotype
(VAF = variant allele frequency)
This question follows from this one.
What is the expected range of values for VAF of mutations with heterozygous genotype (i.e., 0/1 or 0|1)?
Is it possible to have 0% and 100% (respectively having allele-specific depth of 0 in ALT or REF, with non-zero depth in the other allele)?

From the link that I edited into your question (realized that the answer was there):

VAF is the percentage of sequence reads observed matching a specific DNA variant divided by the overall coverage at that locus. Because NGS provides a near random sample, VAF is thus a surrogate measure of the proportion of DNA molecules in the original specimen carrying the variant. For constitutional genetic testing, VAF is a measure of diploid zygosity: heterozygous loci should be near 50% VAF, homozygous loci should be near 100%, and reference loci should be near 0%. Deviations from these three expected values should be considered suspicious as potential errors due to incorrect base calls or alignment. For somatic testing, contamination from normal cells and tumor heterogeneity combine to cause unpredictable VAFs. In some clinical scenarios, such as testing a patient for therapeutic resistance mutations, the desired sensitivity dictates that variants with low VAF be included in downstream analysis. This all means that in order to maintain an acceptable level of sensitivity, variant filtering for somatic variant must be highly permissive from a technical validity perspective.

This is assuming diploid humans. For other ploidies the expected proportions of heterozygosity will match the possible heterozygous genotypes (e.g. 33%/67% for triploid heterozygotes, 25%/75% for tetraploid, etc.).


text:
Why does gene count increase drastically after scaffolding with Hi-C data?
I have a conceptual question.
I have a diploid, outcrossing plant genome assembly of ~1.1 Gb size. The original assembly is generated from PacBio reads. After genome annotation with the Mercator software, the original assembly returned a gene count of 59,615. Then, the diploid genome assembly was phased with Hifiasm+Hi-C data. There, after running the annotation on two haplotype assemblies, the gene count was 75,936 for haplotype 1 and 67,233 for haplotype 2.
Then, we sent the pacbio assembly together with the haplotypes to Dovetail genomics for scaffolding. After the scaffolding with HiRise and annotation with Mercator, I got the following gene count for the three assemblies:
Scaffolded_assembly- 75,874
Scaffolded_haplotype_1- 67,243
Scaffolded_haplotype_2- 86,465

Can anyone please shed some light on why there is such an improvement in gene count? Also, the gene count for haplotypes seems to be the opposite after scaffolding. Do you have any ideas? I do not know how many of those are high confidence or low confidence gene or splice variants.
The BUSCO score was similar in all the cases.
It would be nice if you can give some feedback or any paper explaining such phenomena. Thanks.


text:
How to combine multiple .fasta files of primary assembly from Ensembl into one for sequence alignment?
I have some marmoset snRNA reads that I want to align with the reference transcriptome using cellranger. The primary assembly for marmoset is available here, which is broken down into 22 parts. However, cellranger mkref only accepts one .fa file to generate the transcriptome.
I tried concatentaing all the extracted .fa files into one (After removing the first line which is a comment), but I get the following error while trying to make a reference transcriptome using cellranger:
Creating new reference folder at /data/users/abhishek/ensembl/common_marmoset_ref/combined/common_marmoset
 ...done

 Writing genome FASTA file into reference folder...
 ...done

 Indexing genome FASTA file...
 Traceback (most recent call last):
   File "/data/public-data/cellranger-6.1.1/bin/rna/mkref", line 137, in <module>
     main()
   File "/data/public-data/cellranger-6.1.1/bin/rna/mkref", line 128, in main
     reference_builder.build_gex_reference()
   File "/data/public-data/cellranger-6.1.1/lib/python/cellranger/reference.py", line 630, in build_gex_reference
     self.process_fasta()
   File "/data/public-data/cellranger-6.1.1/lib/python/cellranger/reference.py", line 529, in process_fasta
     subprocess.check_call(["samtools", "faidx", self.fasta_path])
   File "/data/public-data/cellranger-6.1.1/external/anaconda/lib/python3.8/subprocess.py", line 364, in check_call
     raise CalledProcessError(retcode, cmd)
 subprocess.CalledProcessError: Command '['samtools', 'faidx', '/data/users/abhishek/ensembl/common_marmoset_ref/combined/common_marmoset/fasta/genome.fa']' returned non-zero exit     status 1.

There are also files which contain only one or two characters. I suspect that these aren't being handled properly with the direct concatenation of all the files. Some examples are below:
==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.1.fa <==
CCCATTATCCCTGCTTCTGCTAATCTTTCAAAATCTCACCACCATAGTC

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.20.fa <==
TTTTT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.21.fa <==
TT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.22.fa <==
AAAATTGTTTCTTGTTT

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.2.fa <==
C

==> Callithrix_jacchus.mCalJac1.pat.X.dna.primary_assembly.3.fa <==
TTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGG

How do I combine all these .fa files to generate the reference transcriptome for downstream analysis?

When you say "After removing the first line which is a comment", do you mean you removed the lines which start with a >? Those are the headers and are required to have a valid fasta file. They contain the name or identifier of the contig/chromosome/transcript. They separate the contigs in the fasta if there are multiple and will be needed for alignment in order to know which contigs your reads mapped to.

Why do not you run "samtools faidx" alone on the fasta file any see if you are getting any useful error message. Looks like your program is calling samtools for indexing. You can tell this from the last line of the error.

And your fasta file does not look correct.


text:
subsetting more than one ccle sample
I'm working with the CCLE dataset and I'm trying to subset the data to just 3 cell lines of interest but not sure how to it.
Currently, I have been subsetting one by one but is there an easier way?
ccle_panc <- ccle_dataset$data[ , ccle_dataset$colAnnotation$SAMPLE_NAME=="KP4", "", ""]

Assuming you're trying to write R code, array references can use boolean values, and in is useful when you want to select from more than one value:
ccle_panc <- ccle_dataset$data[ , ccle_dataset$colAnnotation$SAMPLE_NAME %in% c("KP4", "DC1", "DC2")]

You can find more information about vector and array subsetting with boolean and numerical values here:
https://datacarpentry.org/R-ecology-lesson/01-intro-to-r.html#Subsetting_vectors

Assuming that you have a dataframe/tibble situation; you can also use the tidyverse approach and the comment from @gringer:
ccle_dataset$data %>%
filter(SAMPLE_NAME %in% c("KP4", "DC1", "DC2"))


text:
Extract only exon regions from GFF/GTF file with input bed regions
I have a bed file that contains a list of regions and I would like to extract only the exon regions from a GTF/GFF file that are contained in the given regions.
Is there a simple command or tool that I can use for that?

You could use BEDOPS:
$ gtf2bed < annotations.gtf | grep -wF exon > exons.bed
$ sort-bed myRegions.unsorted.bed > myRegions.bed
$ bedops --element-of 100% exons.bed myRegions.bed > answer.bed

Or as a one-liner, which is even faster:
$ bedops --element-of 100% <(gtf2bed < annotations.gtf | grep -wF exon) <(sort-bed myRegions.unsorted.bed) > answer.bed

The --element-of 100% argument requires that the exon fall entirely within a region in the myRegions set. This is the most stringent requirement. The most relaxed requirement would be one base of overlap: --element-of 1. You can adjust to either extreme, or to values in between, depending on your criteria.

Thanks Alex !! I just additionally used awk '{ if ($0 ~ "transcript_id") print $0; else print $0" transcript_id \"\";"; }' annotations.gtf as it was giving error Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)

Bedtools intersect will do this:

By far, the most common question asked of two sets of genomic features is whether or not any of the features in the two sets “overlap” with one another. This is known as feature intersection. bedtools intersect allows one to screen for overlaps between two sets of genomic features. Moreover, it allows one to have fine control as to how the intersections are reported. bedtools intersect works with both BED/GFF/VCF and BAM files as input.

https://bedtools.readthedocs.io/en/latest/content/tools/intersect.html


text:
How to calculate switch errors on plink formatted phased genotype data?
everyone!
I'm working with haplotype inference algorithms in pedigree genotype data, and I would like to know if there is a tool or script to calculate switch errors on phased pedigree data (Plink formatted)?
I already know that VCFtools and other scripts do exist, but they only work with VCF format.
Note: I tried to convert my phased pedigree data to VCF but I lost phase, so, this approach does not seem to be helpful, that's why I'm looking for a specific tool to detect switches on phased pedigree data.
A switch errors ocurr when heterozygous positions are misassigned relative to the previous heterozygous position.
See visual example (block 2), True haplotypes up, inferred haplotypes bottom:
https://www.researchgate.net/publication/310766263/figure/fig13/AS:431688409980932@1479933945501/Illustration-of-switch-errors-and-block-breakpoints-In-haplotype-block-2-there-are-two.png
In my case, I got the phased pedigree plink file as follows:
Infered Haplotype, individual 1:
1|2 2|1 1|2 1|2 1|2
Real Haplotype, individual 1:
1|2 1|2 2|1 1|2 1|2
We can see 2 switches, between position 1-2 and 3-4.
Thank you very much in advance!
Regards!

Can you explain what switch errors are?

A switch errors occur when heterozygous positions are misassigned relative to the previous heterozygous position.
See visual example (block 2), True haplotypes up, inferred haplotypes bottom:
https://www.researchgate.net/publication/310766263/figure/fig13/AS:431688409980932@1479933945501/Illustration-of-switch-errors-and-block-breakpoints-In-haplotype-block-2-there-are-two.png
I got the phased pedigree plink file as follows:
Inferred Haplotype, individual 1:
1 1|2 2|1 1|2 1|2 1|2
Real Haplotype, individual 1:
1 1|2 1|2 2|1 1|2 1|2
See 2 switches, between pos. 1-2 and 3-4.


text:
Efficient way of mapping UniProt IDs to representative UniRef90 IDs?
I want to map ~500,000 UniProt IDs to their representative UniRef90 IDs.
What is the most efficient way of doing this?
Do we have a table with the mapping?

You can do this directly on UniProt: https://www.uniprot.org/uploadlists/
Just paste or upload your list of UniProt IDs, and select "UniProtKB AC/ID" in the "From" field and "UniParc" in the "To" field

I've also written a script, pasted below, that can do this with some useful options:
$ uniprot_map.pl -h

uniprot_map.pl will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       uniprot_map.pl [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.

If you save the script as uniprot_map.pl in a directory in your PATH and make it executable, then you can use a file like this:
$ cat file
P31946
P62258

And feed it to the script with the appropriate options to get the mapping:
$ uniprot_map.pl -f ID -t UniParc file
P31946  UPI000013C714
P62258  UPI0000021A46
2 of 2

The script:
#!/usr/bin/perl
use strict;
use warnings;
use LWP::UserAgent;
use Getopt::Std;
my %opts;
getopts('Lsvuhf:F:t:l', \%opts);
my $base = 'http://www.uniprot.org';
my $verbose=$opts{v}||undef;
my $tool = 'mapping';
my (@NAMES,@names);
my $print_only_first=$opts{l}||undef;
my $allow_trembl=$opts{u}||undef;
list_ids() if $opts{L};
&usage() if $opts{h};

&usage() unless $ARGV[0];

####################################
# Collect the IDs we want to match #
####################################
foreach my $ar(@ARGV){
if (-e $ar){
open(A, $ar);
while(<A>){
    chomp;
    push @names,$_;
}
}
else{push @names,$ar};
}

######################################
# What are we matching FROM and TO ? #
######################################
my $from=$opts{f}||'ID';
my $to = $opts{t}||'ACC';

if($opts{t}){
$to=translate_terms($to);
}
if($opts{f}){
    $from=translate_terms($from);
}
#################
# Output format #
#################
my $format = $opts{F}||'tab';
#my @query=

my $agent = LWP::UserAgent->new;
push @{$agent->requests_redirectable}, 'POST';

my $params= {
'from' => $from,
'to' => $to,
'format' => $format,
    'query' => "@names"
};

print STDERR "Requesting...\t" if $verbose;
my $try_again=1;
my $response;
while($try_again==1){
    $response = get_response($params);
    if(defined($response->is_success)){$try_again=0}
elsif($verbose){print STDERR "Failed, retrying...\n"; }
}

if( $print_only_first)
{
if($response->is_success){
my %k;
my @a=split(/\n/,$response->content);
foreach my $line(@a){
    next if $line=~/From/;
    my @b=split(/\t/,$line);
    next if(defined($k{$b[0]}));
    $k{$b[0]}=$b[1];
}
map{print STDOUT "$_\t$k{$_}\n"}keys(%k);
}
else{
die 'Failed, got ' . $response->status_line . ' for ' . $response->request->uri . "\n";
}
}
else{
if($response->is_success){
my %map;
my @a=split(/\n/,$response->content);
#############################
    # remove "From   To"   line #
    #############################
shift(@a) if $a[0] =~/^From/i;
my %kk;
map{
    ########################################################
        # UniProt sometimes returns the SAME line twice    #
        ########################################################
    $kk{$_}++;
    #print STDOUT "$_\n" unless $kk{$_}>1;
    next if $kk{$_}>1;
    my @b=split(/\t/);
    push @{$map{$b[0]}},$b[1];
    }@a;

    foreach my $query (keys(%map)) {
    my %sw;
    foreach my $target (@{$map{$query}}) {
    #############################################
            # Does this prot have SWISSPROT IDs?        #
            #############################################
    if ($to eq 'ID')
    {
        my @b=split(/_/,$target);
        $sw{$target}++ if length($b[0])!=6;
    }
    }
    foreach my $target (@{$map{$query}}) {
    ################################################
    # If a protein has a least one SWISSPROT ID    #
    # skip printing its trEMBL ones unless -u      #
    ################################################
    if ($to eq 'ID' and not $allow_trembl and scalar(keys(%sw)) > 0) {
        print "$query\t$target\n" if defined($sw{$target});
    }
    else {
        print "$query\t$target\n";
        }
        }
    }

    ## Count how many were actually found
    my %found;
    foreach my $resp (@a){
    #print STDERR "R:$resp\t";
    foreach my $id (@names){
    #print STDERR "I:$id\n";
    $found{$id}++ if $resp=~/^$id\t/;
        }

    }
    print STDERR scalar(keys(%found)) . " of " . scalar @names . "\n";

    }
    else{
        die 'Failed, got ' . $response->status_line .
    ' for ' . $response->request->uri . "\n";
    }
}
################################################
# This will translate $to and $from values     #
# between the names used in uniprot flat files #
# and those used by the uniprot database       #
################################################
sub translate_terms{
    my $id=shift();
my %hash=(
"AC" => "ACC", "ID" => "ID", "UniParc" => "UPARC", "UniRef50" => "NF50", "UniRef90" => "NF90", "UniRef100" => "NF100", "EMBL/GenBank/DDBJ" => "EMBL_ID", "EMBL/GenBank/DDBJ CDS" => "EMBLD", "PIR" => "PIR", "UniGene" => "UNIGENE_ID", "GI number*" => "P_GI", "IPI" => "P_IPI", "RefSeq" => "P_REFSEQ_AC", "PDB" => "PDB_ID", "DisProt" => "DISPROT_ID", "HSSP" => "HSSP_ID", "DIP" => "DIP_ID", "MINT" => "MINT_ID", "MEROPS" => "MEROPS_ID", "PeroxiBase" => "PEROXIBASE_ID", "PptaseDB" => "PPTASEDB_ID", "REBASE" => "REBASE_ID", "TCDB" => "TCDB_ID", "Aarhus/Ghent-2DPAGE" => "AARHUS_GHENT_2DPAGE_ID", "ECO2DBASE" => "ECO2DBASE_ID", "World-2DPAGE" => "WORLD_2DPAGE_ID", "Ensembl" => "ENSEMBL_ID", "Ensembl Protein" => "ENSEMBL_PRO_ID", "Ensembl  Trscrpt" => "ENSEMBL_TRS_ID", "Ensembl  Genomes" => "ENSEMBLGENOME_ID", "Ensembl    Genomes" => "Protein", "Ensembl Genomes" => "Transcript", "GeneID" => "P_ENTREZGENEID", "GenomeReviews" => "GENOMEREVIEWS_ID", "KEGG" => "KEGG_ID", "TIGR" => "TIGR_ID", "UCSC" => "UCSC_ID", "VectorBase" => "VECTORBASE_ID", "AGD" => "AGD_ID", "ArachnoServer" => "ARACHNOSERVER_ID", "CGD" => "CGD", "ConoServer" => "CONOSERVER_ID", "CYGD" => "CYGD_ID", "dictyBase" => "DICTYBASE_ID", "EchoBASE" => "ECHOBASE_ID", "EcoGene" => "ECOGENE_ID", "euHCVdb" => "EUHCVDB_ID", "EuPathDB" => "EUPATHDB_ID", "FlyBase" => "FLYBASE_ID", "GeneCards" => "GENECARDS_ID", "GeneDB_Spombe" => "GENEDB_SPOMBE_ID", "GeneFarm" => "GENEFARM_ID", "GenoList" => "GENOLIST_ID", "H-InvDB" => "H_INVDB_ID", "HGNC" => "HGNC_ID", "HPA" => "HPA_ID", "LegioList" => "LEGIOLIST_ID", "Leproma" => "LEPROMA_ID", "MaizeGDB" => "MAIZEGDB_ID", "MIM" => "MIM_ID", "MGI" => "MGI_ID", "NMPDR" => "NMPDR", "Orphanet" => "ORPHANET_ID", "PharmGKB" => "PHARMGKB_ID", "PseudoCAP" => "PSEUDOCAP_ID", "RGD" => "RGD_ID", "SGD" => "SGD_ID", "TAIR" => "TAIR_ID", "TubercuList" => "TUBERCULIST_ID", "WormBase" => "WORMBASE_ID", "WormBase Trscrpt" => "WORMBASE_TRS_ID", "WormBase Protein" => "WORMBASE_PRO_ID", "Xenbase" => "XENBASE_ID", "ZFIN" => "ZFIN_ID", "eggNOG" => "EGGNOG_ID", "HOGENOM" => "HOGENOM_ID", "HOVERGEN" => "HOVERGEN_ID", "OMA" => "OMA_ID", "OrthoDB" => "ORTHODB_ID", "ProtClustDB" => "PROTCLUSTDB_ID", "BioCyc" => "BIOCYC_ID", "Reactome" => "REACTOME_ID", "CleanEx" => "CLEANEX_ID", "GermOnline" => "GERMONLINE_ID", "DrugBank" => "DRUGBANK_ID", "NextBio" => "NEXTBIO_ID"
);
defined($hash{$id}) ?
return($hash{$id}) :
return($id) ;
}

sub get_response{
    my $params=shift;
my $response = $agent->post("$base/$tool/", $params);
    while (my $wait = $response->header('Retry-After')) {
    print STDERR "Waiting ($wait)...\n";
sleep $wait;
    $response = $agent->get($response->base);
}
return($response);

}

sub usage{
    $0=~s/.+?([^\/]+)$/$1/;

    print STDERR <<EndOfHelp;

$0 will connect to UniProt and map IDs from the input file to ehatever ID was requested.

USAGE:
       $0 [options] <ID_LIST>

OPTIONS:
      -h : Print this help and exit.
      -f : From ID type. This is the type of ID found in the input file.
      -F : Format. See http://www.uniprot.org/faq/28#batch_retrieval_of_entries for details.

      -l : When a requested ID maps to many output IDs, print only the first (probably the best) match.
      -L : List all recognised ID types.
      -u : Print unreviewed (trEMBL) prots as well. By default, if an ID maps to SWISSPROT and trEMBL IDs
           only the SWISSPROTs will be printed, use this flag to print trEMBL as well.
      -t : To ID. This is the ID type that the input IDs will be mapped to.
      -v : Verbose. Print progress messages to STDERR.
EndOfHelp

    exit(0);

}
sub list_ids{
    open(HELP, "| less") ;
    print HELP <<EndOfHelp;

Database            Abbreviation
===================== UniProt =============================
UniProtKB AC/ID :   ACC+ID
UniProtKB AC    :   ACC
UniProtKB ID    :   ID
UniParc         :   UPARC
UniRef50    :   NF50
UniRef90    :   NF90
UniRef100   :   NF100

=============== Other sequence databases ==================
EMBL/GenBank/DDBJ:  EMBL_ID
EMBL/GenBank/DDBJ CDS:  EMBLD
PIR     :   PIR
UniGene     :   UNIGENE_ID
GI number*  :   P_GI
IPI     :   P_IPI
RefSeq      :   P_REFSEQ_AC
HGNC name       :       GENENAME
================== 3D structure databases =================
PDB     :   PDB_ID
DisProt     :   DISPROT_ID
HSSP        :   HSSP_ID

============ Protein-protein interaction databases ========
DIP     :   DIP_ID
MINT        :   MINT_ID

============ Protein family/group databases ===============
MEROPS      :   MEROPS_ID
PeroxiBase  :   PEROXIBASE_ID
PptaseDB    :   PPTASEDB_ID
REBASE      :   REBASE_ID
TCDB        :   TCDB_ID

================== 2D gel databases =======================
Aarhus/Ghent-2DPAGE:  AARHUS_GHENT_2DPAGE_ID
ECO2DBASE   :   ECO2DBASE_ID
World-2DPAGE    :   WORLD_2DPAGE_ID

============= Genome annotation databases =================
Ensembl         :   ENSEMBL_ID
Ensembl Protein :   ENSEMBL_PRO_ID
Ensembl Trscrpt :   ENSEMBL_TRS_ID
Ensembl Genomes :   ENSEMBLGENOME_ID
Ensembl Genomes :   Protein
Ensembl Genomes :   Transcript
GeneID          :   P_ENTREZGENEID
GenomeReviews   :   GENOMEREVIEWS_ID
KEGG        :   KEGG_ID
TIGR        :   TIGR_ID
UCSC        :   UCSC_ID
VectorBase  :   VECTORBASE_ID

============ Organism-specific gene databases =============
AGD     :   AGD_ID
ArachnoServer   :   ARACHNOSERVER_ID
CGD     :   CGD
ConoServer  :   CONOSERVER_ID
CYGD        :   CYGD_ID
dictyBase   :   DICTYBASE_ID
EchoBASE    :   ECHOBASE_ID
EcoGene     :   ECOGENE_ID
euHCVdb     :   EUHCVDB_ID
EuPathDB    :   EUPATHDB_ID
FlyBase     :   FLYBASE_ID
GeneCards   :   GENECARDS_ID
GeneDB_Spombe   :   GENEDB_SPOMBE_ID
GeneFarm    :   GENEFARM_ID
GenoList    :   GENOLIST_ID
H-InvDB     :   H_INVDB_ID
HGNC        :   HGNC_ID
HPA     :   HPA_ID
LegioList   :   LEGIOLIST_ID
Leproma     :   LEPROMA_ID
MaizeGDB    :   MAIZEGDB_ID
MIM     :   MIM_ID
MGI     :   MGI_ID
NMPDR       :   NMPDR
Orphanet    :   ORPHANET_ID
PharmGKB    :   PHARMGKB_ID
PseudoCAP   :   PSEUDOCAP_ID
RGD     :   RGD_ID
SGD     :   SGD_ID
TAIR        :   TAIR_ID
TubercuList :   TUBERCULIST_ID
WormBase    :   WORMBASE_ID
WormBase Trscrpt:   WORMBASE_TRS_ID
WormBase Protein:   WORMBASE_PRO_ID
Xenbase     :   XENBASE_ID
ZFIN        :   ZFIN_ID

================= Phylogenomic databases ==================
eggNOG      :   EGGNOG_ID
HOGENOM     :   HOGENOM_ID
HOVERGEN    :   HOVERGEN_ID
OMA     :   OMA_ID
OrthoDB     :   ORTHODB_ID
ProtClustDB :   PROTCLUSTDB_ID

=================== Enzyme and pathway ====================
BioCyc      :   BIOCYC_ID
Reactome    :   REACTOME_ID

================= Gene expression databases ===============
CleanEx     :   CLEANEX_ID
GermOnline  :   GERMONLINE_ID

========================== Other ==========================
DrugBank    :   DRUGBANK_ID
NextBio     :   NEXTBIO_ID

EndOfHelp
close(HELP);
    exit(0);
}

If you have too many IDs to do in a single batch, you could either simply split the file into smaller ones and run the script on each file in a loop or modify the script to do it.
Split the input:
$ wc -l names
$ wc -l names
500000 names
## split into 500 files of 1K lines each
$ split -l 1000 names chunk
## iterate over each chunk of names
$ for i in chunk*; do uniprot_map.pl -f ID -t UniParc "$i" >> mapped; done

Thank you @terdon for the detailed answer and your script! Let me try it out!

Seems both of the solutions require ID list < 2Mb. I will try to look into the XML file of UniProt. Please check out: https://blog.liang2.tw/posts/2018/01/read-uniprotkb-xml/

@JohnnyTam just split your file and run the script in a loop. Or change the script to run them in batches.

@JohnnyTam see updated answer for how to split the input. Alternatively, if you can give me your list of names, I can have a go at modifying the script to handle them.

Ops, I forgot I can split. It worked perfectly! Thank you very much!


text:
Why does this modification site file contain duplicate sites and multiple related proteins?
The data below is coming from a MaxQuant post-translational protein modification (PTM) site file.
Why does it contain duplicate sites (yellow) and proteins (cyan)?

Source data
http://ftp.ebi.ac.uk/pride-archive/2020/12/PXD020517/MaxQuant_BA39_SOLUBLE.zip

There are multiple protein isoforms that the gene may encode. The hyphen-number after the Uniprot ID represents the isoform —in a Uniprot you find the detail of these at the bottom.
The canonical isoform in Uniprot is the most likely one to be biologically relevant, while in NCBI and Ensembl it is the longest.
In a limited amount of case multiple isoforms are biologically relevant.
For the sake of sanity most variant annotation tools have a canonical argument to use the canonical.

According to a proteomics researcher at the NIH:

proteins column contains uniprot IDs.
position within proteins column are the starting position of each ID.
position is the leading protein


text:
Under what circumstances can the SEQ field in a SAM file be a *?
I saw in the SAM format specifications that the SEQ field (10th column) can be a "*" if the sequence is not stored, instead of being the sequence of the mapped read. Under what circumstances is this expected to occur? I don't understand why it's allowed in the spec.
For a concrete example, I have pacbio data mapped to a reference sequence, and some of the records in the SAM file have a sequence of *, and I don't know why or if that means anything. I don't understand why they exist, since the data is ostensibly present in the FASTQ file I use as input to the read mapper, and I don't understand what purpose they could have, since it seems to me that these records couldn't be interpreted by variant calling software. At best it would seem they can be used to infer coverage, but it's coverage by sequences of an uncertain nature?

Related: How to output all sequences with bwa mem, not *?

Sounds like I should check if these are exclusively secondary alignments

This could be a duplicate question as @acvill is indicating


text:
What are values of FILTER column of vcf files produced by Mutect2
I have called SNVs with Mutect2, now in filter column of a file called filtered vcf I have a lot of things like
> length(unique(filtered.vcf$FILTER)) [1] 409 >

> unique(filtered.vcf $FILTER)
  [1] "germline"
  [2] "PASS"
  [3] "clustered_events"
  [4] "clustered_events;panel_of_normals"
  [5] "weak_evidence"
  [6] "haplotype"
  [7] "base_qual;haplotype"
  [8] "clustered_events;germline;haplotype"
  [9] "base_qual;clustered_events;germline;haplotype"
 [10] "clustered_events;haplotype"
 [11] "clustered_events;haplotype;weak_evidence"
 [12] "base_qual;clustered_events;haplotype"
 [13] "haplotype;map_qual"
 [14] "germline;panel_of_normals"
 [15] "base_qual"
 [16] "germline;haplotype"
 [17] "germline;panel_of_normals;position"
 [18] "clustered_events;orientation"
 [19] "clustered_events;orientation;strand_bias"
 [20] "base_qual;haplotype;map_qual"
 [21] "clustered_events;fragment;haplotype"
 [22] "base_qual;clustered_events"
 [23] "germline;haplotype;panel_of_normals"
 [24] "clustered_events;germline;haplotype;panel_of_normals"
 [25] "fragment"
 [26] "clustered_events;haplotype;map_qual"
 [27] "clustered_events;germline"
 [28] "clustered_events;germline;panel_of_normals"
 [29] "germline;weak_evidence"
 [30] "orientation"

I tought that the filtered.vcf Mutect2 output file contain only the filtered events, but the values I see suggest otherwise. How should I understand the output?

You haven't told us what commands you ran, so I am assuming you first ran Mutect2 and then FilterMutectCalls. If so, the FilterMutectCalls does not remove any variants from its input, it will simply flag any that fail its builtin filters. Any variants that did not fail a filter are marked as PASS and the rest have a filter added.
Now, whether you want to filter these out or not is up to you. You need to look at the VCF header where each of the filters is defined (admittedly, not very clearly) and make a judgment call on whether you consider this filter to be a deal breaker or not.
For instance, the germline flag means:
##FILTER=<ID=germline,Description="Evidence indicates this site is germline, not somatic">

So that indicates that gatk believes this is more likely to be a germline variant than a somatic one, that it is more likely to be a constituent variant present in the sample's germline cells than a variant that is only found in the sample's tumor cells.
Another clear case is the base_qual filter which means:
##FILTER=<ID=base_qual,Description="alt median base quality">

When a variant is flagged with this filter, it means that the median fastq base quality of the reads supporting the variant was significantly lower than the median fastq base quality of the reads supporting the reference alele at the same position, making it likely that the variant is a false positive.
You need to read the VCF headers and any gatk documentation you can find (warning: these filters are not very well documented at all, in my experience), understand what the filters are and then decide what variants you consider real based on what you know about your sample, your experimental design and the question you are trying to answer.

How did you get to the explanation of base_qual from just the "alt median base quality"? I'd like to explain each filter in simple English like you've done for base_qual above and there seems to be very little proper documentation.

@RamRS yes, extremely little. I got this one from parsing  "alt median base quality" as "median base quality in the reads supporting the alternative allele" and it has also been borne out by manual testing. Unfortunately, there is almost no documentation available at all for these things.


text:
How to Sort and Index a SAM file without converting it to BAM?
Generally, I use samtools to sort and index BAM files. Samtools works also with BGZF compressed SAM but not with plain ones.
I have terabytes of plain SAMs at the moment and I would love to avoid compression/conversion. I have all the data in several external HD and I am not allowed to delete the original SAM files. I would need much more space than I have now to also store the BAM files.
So, I was wondering if there is a way to avoid it. Is there some tool/trick to sort and index uncompressed SAM files, without converting them to BAM or compressing them?

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.

In general, avoid storing any SAM files on disk. Instead, use htslib/pysam/etc to directly read BAM files, or pipe SAM to other tools or your own scripts.

Not only will you save disk space by converting to BAM, but BAM files are faster to manipulate than SAM. Source: Dave Tang's SAMTools wiki.
sort supports uncompressed SAM format from a file or stdin, though index requires BGZIP-compressed SAM or BAM. I don't think you can get around this.
If you're dead-set on reading from and writing to uncompressed SAM format, you can use pipes to avoid writing BAM files. Specify the output of view to be an uncompressed bam, then pipe it to your next command, then pipe it to view again to output your sorted SAM file.
samtools view -u in.sam | samtools sort | samtools view -h > sort.sam
samtools view -u sort.sam | samtools index

I suspect you could convince tabix to index a SAM file, though you're correct that using BAM is a MUCH better solution.

Tabix has a preset for sam -p sam, but yeah, use bam, there is no advantage using sam anyway and barely any downstream tool swallows sam.gz.