Convert bam file to highly compressible bam I have a large collection of bam files and I want to post-process each of them into another bam where I can make queries about:

the reads position and pair-endness,
insert sizes,
MAPQ and other flags, etc. of the reads

but where I don't need to make use of:

the sequence string of the reads,
the Q-scores qualities or,
the read ids, apart from keeping pair-endness in the read id,

So I would like to know what's a good tool/parameters to process each bam and turn them into a small size bam that still retains the information needed above. CRAM format is not an option at this point as several of the tools I use are not CRAM-compatible yet.

Any ideas?

You can just set the fields you don't need to *:

samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam -

This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields.

The resulting BAM file looks like:

*   83  chr1    3003308 37  51M =   3003157 -202    *   *
*   99  chr1    3003364 33  51M =   3003547 234 *   *
*   99  chr1    3003410 37  50M =   3003636 269 *   *
*   99  chr1    3003438 255 50M =   3003760 373 *   *


Sibling vs Parent/Child relationship detection using whole-genome sequence data I am analyzing a family-based whole-genome sequence (WGS) dataset.  My identity-by-descent(IBD) analysis shows they share about 50% of their DNA. I think my pedigree structure may have gotten mixed up though and I want to see if one individual is the mother or the sister of a male individual.  Is there anyway to do this either with Plink or KING?

With at least three samples, you can determine this by looking at the IBS0 output column (from either plink2 --make-king-table or KING): if the relationship in question has a very low IBS0 value compared to the other pairs, that's parent-child.  You can try to do the same thing with only two samples, but with no other IBS0 values to compare against it's easier to misjudge.


Why is bwa-mem the standard algorithm when using bwa? The industry standard for aligning short reads seems to be bwa-mem. However, in my tests I have seen that using bwa backtrack (bwa-aln + bwa-sampe + bwa-samse) performs better. It is slightly slower, but gives significantly better results in terms of both sensitivity and specificity. I have tested it using the genome in a bottle data and public samples (NA12878 and NA12877 among others) and found that backtrack consistently outperformed bwa-mem.

So why is bwa-mem the standard? Am I wrong and other tests have shown the opposite? I don't really see how since I tested using the most common datasets and validation data. Is it that the slight increase in efficiency outweighs the decrease in performance?

The only other explanation I can see is that bwa backtrack is designed specifically for Illumina reads and all my tests have been on Illumina data. Is it just that bwa-mem is "sequencer agnostic"? So that we can use the same algorithm irrespective of what sequencing platform is used? In that case, it makes sense to use backtrack for if we only deal with Illumina data and mem if we can have different sequencers. But, if so, seeing as Illumina is so widespread, why isn't backtrack used more often on Illumina data? I feel I must be missing something.

bwa mem is newer, faster, and [should be] more accurate, particularly for longer reads.

From the bwa man page (presumably in Heng Li's own words):

  BWA is a software package for mapping low-divergent sequences against a large reference genome, such as the human genome. It consists of three algorithms:  BWA-backtrack, BWA-SW and BWA-MEM. The first algorithm  is designed for Illumina sequence reads up to 100bp, while the rest two for longer sequences ranged from 70bp to 1Mbp. BWA-MEM and BWA-SW share similar  features such as long-read support and split alignment, but BWA-MEM, which is the latest, is generally recommended for high-quality queries as it is faster and  more  accurate. BWA-MEM also has better performance than BWA-backtrack for 70-100bp Illumina reads.


How do you map nanopore fast5 files? Is there a best practices for mapping Oxford Nanopore files to a reference?

Is there a tool that can take a tarball of fast5 files and map them directly or do they need to be converted to fastq first?

Fast5 files need to be converted to fastq first prior to mapping. There are currently no tools available to map the entire raw signal directly to a reference, although a few tools are available for identifying short sequences within the raw signal.


How to select a cutoff for interaction confidence in STRINGdb? I have a list of 100 genes that are called as hits in a genetic screening. I want to have a network of the interactions between the proteins of these 100 genes. I am using both STRINGdb web and its R API.

In both situations you can select the threshold of confidence of interaction (default 0.4) for displaying the network. Changing this parameter drastically changes the network that I am generating. I have read the FAQ section of STRINGdb and they recommend to choose some arbitrary number based on the number of interactions you need for you analysis. If I use the default threshold, I have few interactions and I don't know if lowering the threshold would be correct.

Therefore, my question is, is there any established threshold based on habitual user's experience used for that, beyond the recommendation of setting an arbitrary one?

This depends on what you are trying to do and whether you value specificity over sensitivity. We can't tell you since it is entirely dependent on the biological question you want to answer.
However, I would recommend two things:

Don't use STRING. The creators of STRING made the choice to value sensitivity over all else, so they include any interaction they can get their hands on. That is sometimes fine, depending on what you want to do, but is more often a problem. Personally, I tend to avoid STRING as much as possible and instead use more curated databases like APID or IntAct. You can find a list of many PPI databases in the EBI's Psiqcuick View page. You might also be interested in my answer here which gives an example script using APID which can easily be modified to query other DBs.

Whatever you use, I recommend you filter by detection method. There are many different interaction detection methods ranging from those that only identify direct, binary interactions between two proteins (e.g. Yeast two-hybrid), those that also find proteins that are in the same complex even if there is no direct interaction between them (e.g. ChIP), to various non-experimental methods which are used to infer interactions.
Luckily, all of these have been collected into a controlled vocabulary which you can see in the Interaction Detection Method Ontology Lookup Service page. Go there, collect the interaction detection methods you consider good enough for whatever you want to do and filter your interactions using those.

For example, for a high quality human interactome I built for a project I used to work on, I was only interested in direct binary interactions, so I only kept interactions detected by the following methods:
    "MI:0008" => "array technology",
    "MI:0009" => "bacterial display",
    "MI:0010" => "beta galactosidase complementation",
    "MI:0011" => "beta lactamase complementation",
    "MI:0012" => "bioluminescence resonance energy transfer",
    "MI:0013" => "biophysical",
    "MI:0014" => "adenylate cyclase complementation",
    "MI:0016" => "circular dichroism",
    "MI:0017" => "classical fluorescence spectroscopy",
    "MI:0018" => "two hybrid",
    "MI:0020" => "transmission electron microscopy",
    "MI:0030" => "cross-linking study",
    "MI:0031" => "protein cross-linking with a bifunctional reagent",
    "MI:0034" => "display technology",
    "MI:0040" => "electron microscopy",
    "MI:0041" => "electron nuclear double resonance",
    "MI:0042" => "electron paramagnetic resonance",
    "MI:0043" => "electron resonance",
    "MI:0047" => "far western blotting",
    "MI:0048" => "filamentous phage display",
    "MI:0049" => "filter binding",
    "MI:0051" => "fluorescence technology",
    "MI:0052" => "fluorescence correlation spectroscopy",
    "MI:0053" => "fluorescence polarization spectroscopy",
    "MI:0055" => "fluorescent resonance energy transfer",
    "MI:0065" => "isothermal titration calorimetry",
    "MI:0066" => "lambda phage display",
    "MI:0073" => "mrna display",
    "MI:0081" => "peptide array",
    "MI:0084" => "phage display",
    "MI:0089" => "protein array",
    "MI:0090" => "protein complementation assay",
    "MI:0091" => "chromatography technology",
    "MI:0092" => "protein in situ array",
    "MI:0095" => "proteinchip(r) on a surface-enhanced laser desorption/ionization",
    "MI:0097" => "reverse ras recruitment system",
    "MI:0098" => "ribosome display",
    "MI:0099" => "scintillation proximity assay",
    "MI:0107" => "surface plasmon resonance",
    "MI:0108" => "t7 phage display",
    "MI:0111" => "dihydrofolate reductase reconstruction",
    "MI:0112" => "ubiquitin reconstruction",
    "MI:0114" => "x-ray crystallography",
    "MI:0115" => "yeast display",
    "MI:0226" => "ion exchange chromatography",
    "MI:0227" => "reverse phase chromatography",
    "MI:0231" => "mammalian protein protein interaction trap",
    "MI:0232" => "transcriptional complementation assay",
    "MI:0255" => "post transcriptional interference",
    "MI:0369" => "lex-a dimerization assay",
    "MI:0370" => "tox-r dimerization assay",
    "MI:0397" => "two hybrid array",
    "MI:0398" => "two hybrid pooling approach",
    "MI:0399" => "two hybrid fragment pooling approach",
    "MI:0400" => "affinity technology",
    "MI:0401" => "biochemical",
    "MI:0405" => "competition binding",
    "MI:0406" => "deacetylase assay",
    "MI:0410" => "electron tomography",
    "MI:0411" => "enzyme linked immunosorbent assay",
    "MI:0415" => "enzymatic study",
    "MI:0416" => "fluorescence microscopy",
    "MI:0419" => "gtpase assay",
    "MI:0420" => "kinase homogeneous time resolved fluorescence",
    "MI:0423" => "in-gel kinase assay",
    "MI:0424" => "protein kinase assay",
    "MI:0425" => "kinase scintillation proximity assay",
    "MI:0426" => "light microscopy",
    "MI:0428" => "imaging technique",
    "MI:0432" => "one hybrid",
    "MI:0434" => "phosphatase assay",
    "MI:0435" => "protease assay",
    "MI:0437" => "protein three hybrid",
    "MI:0440" => "saturation binding",
    "MI:0508" => "deacetylase radiometric assay",
    "MI:0509" => "phosphatase homogeneous time resolved fluorescence",
    "MI:0510" => "homogeneous time resolved fluorescence",
    "MI:0511" => "protease homogeneous time resolved fluorescence",
    "MI:0512" => "zymography",
    "MI:0513" => "collagen film assay",
    "MI:0514" => "in gel phosphatase assay",
    "MI:0515" => "methyltransferase assay",
    "MI:0516" => "methyltransferase radiometric assay",
    "MI:0655" => "lambda repressor two hybrid",
    "MI:0657" => "systematic evolution of ligands by exponential enrichment",
    "MI:0678" => "antibody array",
    "MI:0695" => "sandwich immunoassay",
    "MI:0696" => "polymerase assay",
    "MI:0726" => "reverse two hybrid",
    "MI:0727" => "lexa b52 complementation",
    "MI:0728" => "gal4 vp16 complementation",
    "MI:0809" => "bimolecular fluorescence complementation",
    "MI:0813" => "proximity enzyme linked immunosorbent assay",
    "MI:0824" => "x-ray powder diffraction",
    "MI:0825" => "x-ray fiber diffraction",
    "MI:0827" => "x-ray tomography",
    "MI:0841" => "phosphotransferase assay",
    "MI:0870" => "demethylase assay",
    "MI:0872" => "atomic force microscopy",
    "MI:0879" => "nucleoside triphosphatase assay",
    "MI:0880" => "atpase assay",
    "MI:0887" => "histone acetylase assay",
    "MI:0889" => "acetylase assay",
    "MI:0892" => "solid phase assay",
    "MI:0894" => "electron diffraction",
    "MI:0895" => "protein kinase A complementation",
    "MI:0899" => "p3 filamentous phage display",
    "MI:0900" => "p8 filamentous phage display",
    "MI:0905" => "amplified luminescent proximity homogeneous assay",
    "MI:0916" => "lexa vp16 complementation",
    "MI:0920" => "ribonuclease assay",
    "MI:0921" => "surface plasmon resonance array",
    "MI:0946" => "ping",
    "MI:0947" => "bead aggregation assay",
    "MI:0949" => "gdp/gtp exchange assay",
    "MI:0953" => "polymerization",
    "MI:0968" => "biosensor",
    "MI:0969" => "bio-layer interferometry",
    "MI:0972" => "phosphopantetheinylase assay",
    "MI:0976" => "total internal reflection fluorescence spectroscopy",
    "MI:0979" => "oxidoreductase assay",
    "MI:0984" => "deaminase assay",
    "MI:0989" => "amidase assay",
    "MI:0990" => "cleavage assay",
    "MI:0991" => "lipid cleavage assay",
    "MI:0992" => "defarnesylase assay",
    "MI:0993" => "degeranylase assay",
    "MI:0994" => "demyristoylase assay",
    "MI:0995" => "depalmitoylase assay",
    "MI:0996" => "deformylase assay",
    "MI:0997" => "ubiquitinase assay",
    "MI:0998" => "deubiquitinase assay",
    "MI:0999" => "formylase assay",
    "MI:1000" => "hydroxylase assay",
    "MI:1001" => "lipidase assay",
    "MI:1002" => "myristoylase assay",
    "MI:1003" => "geranylgeranylase assay",
    "MI:1004" => "palmitoylase assay",
    "MI:1005" => "adp ribosylase assay",
    "MI:1006" => "deglycosylase assay",
    "MI:1007" => "glycosylase assay",
    "MI:1008" => "sumoylase assay",
    "MI:1009" => "desumoylase assay",
    "MI:1010" => "neddylase assay",
    "MI:1011" => "deneddylase assay",
    "MI:1016" => "fluorescence recovery after photobleaching",
    "MI:1019" => "protein phosphatase assay",
    "MI:1024" => "scanning electron microscopy",
    "MI:1026" => "diphtamidase assay",
    "MI:1030" => "excimer fluorescence",
    "MI:1031" => "protein folding/unfolding",
    "MI:1036" => "nucleotide exchange assay",
    "MI:1037" => "Split renilla luciferase complementation",
    "MI:1038" => "silicon nanowire field-effect transistor",
    "MI:1087" => "monoclonal antibody blockade",
    "MI:1088" => "phenotype-based detection assay",
    "MI:1089" => "nuclear translocation assay",
    "MI:1111" => "two hybrid bait or prey pooling approach",
    "MI:1112" => "two hybrid prey pooling approach",
    "MI:1113" => "two hybrid bait and prey pooling approach",
    "MI:1137" => "carboxylation assay",
    "MI:1138" => "decarboxylation assay",
    "MI:1142" => "aminoacylation assay",
    "MI:1145" => "phospholipase assay",
    "MI:1147" => "ampylation assay"

If you also want those that might not be direct binary but simply in the same complex, include:
      "MI:0019" => "coimmunoprecipitation",
      "MI:0006" => "anti bait coimmunoprecipitation",
      "MI:0007" => "anti tag coimmunoprecipitation",
      "MI:0858" => "immunodepleted coimmunoprecipitation",
      "MI:0096" => "pull down",
      "MI:0963" => "interactome parallel affinity capture",
      "MI:0676" => "tandem affinity purification"


Is there a difference between any scRNA-seq visualization and pseudotime? Is there a good definition of pseudotime? Some tools are clearly labeled as pseudotime and produce values along a trajectory, but there are more complex approaches that involve branching and disconnected populations. There are 70 of them in the dynverse, so there are many options.
UMAP will also order cells, but that is generally not considered pseudotime. You can even overlay Monocle or PAGA results onto a UMAP so the cells are still ordered the same way. SPRING is specifically designed to capture continuous trajectories, but is also not considered pseudotime. Why not?

Pseudotime is the order of cells along a trajectory which in the simplest case is a minimum-spanning tree based on a (low-dimensional) manifold such as PCA or UMAP. Without a defined trajectory you have no pseudotime, and therefore UMAP is not pseudotime, but in fact can be used to define the trajectory.


How can I translate infix to mathml using the python interface to libsbml? Assume I have an infix expression such as k1*s1, how do I use libsbml to translate this expression to mathml using the python interface?

It’s quite straightforward to convert infix to mathml using libsbml. Make sure you have libsbml installed in python by typing

pip install python-libsbml

at the windows/Mac/Linux command line or

!pip install python-libsbml

if you are running a Jupyter notebook. Some tools such as Tellurium already have libsbml preinstalled.

At the python command line type

import libsbml

to import the package. To enter the infix expression type the following:

p = libsbml.parseFormula ('k1*s1')

To generate the mathml type

print (libsbml.writeMathMLToString (p))


snakemake get environment variable for rule name and log path When a rule fails in snakemake workflow, it returns a descriptive message in log file like following:

Error in rule step1:
    jobid: 99
    output: results/foo/step1/sample5/outdir/output1.tsv, results/foo/step1/sample5/outdir/output2.tsv
    log: logs/step1/sample5_out.log (check log file(s) for error message)
    shell:
    ...
        (exited with non-zero exit code)

Since I run workflow over many samples, I have defined a postrun.sh bash script in jobscript.sh directive that runs after snakemake {exec_job} section and notifies me via userbot in slack channel when a step failed for one or more samples.

jobscript.sh

        ########################### START SNAKEMAKE CMD ############################
    {exec_job}
        ############################ END SNAKEMAKE CMD #############################
    exitstat=$?
    export exitstat

    echo "END at $(date)"
. "${{HOME}}"/.smk_confs/postrun.sh

postrun.sh

## execute after running each of snakemake job
echo 'Running postrun.sh'
ERRMSG="MYJOB ID: $PBS_JOBID failed at $(pwd) on $(hostname) for ${USER} with exit status: ${exitstat}."

## ping to slack
ping_slack -m "$ERRMSG"
sleep 5

I am looking for a snakemake exported environment variable that I can use in postrun.sh to know which rule ID failed for one or more samples. For now, I can think of using job scheduler variable matching job name, e.g, PBS_JOBNAME, but better if I have a variable from snakemake API for rule ID or preferably, path to log file.

Thanks

You could add an exit trap to the beginning of the shell scripts inside your Snakefile, which is executed after each script is done and sends a message when there was an error:

 shell:
    """
    trap 'if [ $? -neq 0 ]; then ping_slack -m "Error in rule xxx with input {input}"; exit 1' EXIT
    ...


What is the fastest way to calculate the number of unknown nucleotides in FASTA / FASTQ files? I used to work with publicly available genomic references, where basic statistics are usually available and if they are not, you have to compute them only once so there is no reason to worry about performance.

Recently I started sequencing project of couple of different species with mid-sized genomes (~Gbp) and during testing of different assembly pipelines I had compute number of unknown nucleotides many times in both raw reads (in fastq) and assembly scaffolds (in fasta), therefore I thought that I would like to optimize the computation.

For me it is reasonable to expect 4-line formatted fastq files, but general solution is still prefered
It would be nice if solution would work on gzipped files as well

Q : What is the fastest way (performance-wise) to compute the number of unknown nucleotides (Ns) in fasta and fastq files?


What is local realignment and what is the problem it solves? I am trying to understand local realignment but I could not get a clear idea of what is the problem solved by it.

For example, reading Homer and Nelson (2010):

  Because alignment algorithms map each read independently to the reference genome, alignment artifacts could result, such that SNPs, insertions, and deletions are improperly placed relative to their true location. This leads to local alignment errors due to a combination of sequencing error, equivalent positions of the variant being equally likely, and adjacent variants or nearby errors driving misalignment of the local sequence. These local misalignments lead to false positive variant detection, especially at apparent heterozygous positions. For example, insertions and deletions towards the ends of reads are difficult to anchor and resolve without the use of multiple reads. In some cases, strict quality and filtering thresholds are used to overcome the false detection of variants, at the cost of reducing power [13]. Since each read represents an independent observation of only one of two possible haplotypes (assuming a diploid genome), multiple read observations could significantly reduce false-positive detection of variants.

I don't understand how the different artifacts (such as SNPs) can result from the alignment.

Homer, Nils, and Stanley F. Nelson. "Improved variant discovery through local re-alignment of short-read next-generation sequencing data using SRMA." Genome biology 11.10 (2010): R99.

SNPs are likely to be created and InDels are likely to be missed. Suppose you have a read, ACTGACTGACTGTAC and you align it to a reference sequence ACTGACTGACTGTTAAGAACGACTACGAC. If you aligned that, you would either get:

ACTGACTGACTGTac    (lower case denotes soft-clipping)
ACTGACTGACTGTTAAGAACGACTACGAC

or

ACTGACTGACTGTAC    (N.B., you've created some possible SNPS)
ACTGACTGACTGTTAAGAACGACTACGAC

But what happens if you look at a bunch of reads? See below:

ACTGACTGACTGTAC
  TGACTGACTGTACTA
    ACTGACTGTACTAAG
ACTGACTGACTGT--TAAGAACGACTACGAC

I've made the 2 base insertion obvious. As you can see, using many alignments, we can "realign" our original read to see that it actually ends in a 2 base insertion, something that's completely impossible to find with a single read. This realignment is part of what the GATK haplotype caller does (it makes a local de Bruijn graph). That's also what my MethIndelRealigner does. I have another real-world example in its documentation.

Raw alignments in IGV:

Alignments after realignment:

One can find plenty of similar scenarios that result in SNPs rather than simply losing coverage for InDels.


Difference between samtools mark duplicates and samtools remove duplicates? What is difference between samtools mark duplicates and remove duplicates ?  Is it necessary to mark duplicates before removing duplicates with samtools?

samtools rmdup and samtools markdup -r do the same thing. Without the -r flag samtools markdup only flags the duplicates.

You'll have to run samtools fixmate -m and sort the output to add ms and MC tags to prepare the file for markdup..

According to the current documentation, rmdup is obsolete. Please do not use rmdup.


Extract end position from CIGAR I was wondering how could I obtain the end position of an alignment using its start position and its CIGAR String.

If an alignment present some soft clipping, for example:

    20S55M

Then, would the start position indicated in the SAM file be the first soft clipped base, or the first correctly mapped base?

Thanks

The start position is the first correctly mapped base. One way to do it, is to use GenomicAlignments in R to read in the alignments.

For example the sam file looks like this:

cat test.sam
@HD VN:1.3  SO:coordinate
@SQ SN:chrI LN:228861
D00689:415:CD4LUANXX:3:1316:7447:98083#CGGCTATGCAGGACGT 0   chrI    769 0   42M8S   *   0   0   CCATAGCCCGGTGGGGCCTGTGCTGCCACAGAGGGTGTGTCTAGATCGGA  B/B//FFFFFFBFFFF//7FF/<F/FFBFFFF//<F///<F/BFB//<F/  NM:i:1  MD:Z:21A20  AS:i:37 XS:i:37
D00689:415:CD4LUANXX:3:2308:9597:30445#CGGCTATGCAGGACGT 0   chrI    3965    0   3S47M   *   0   0   CCACACACCCACACACCACACCCACACCACACCCACACCACACCCACACC  BBBBBFFFFFFFFFFFFFFFFFFFFFFFFFFFFBFFFFFBF<FFFFF<FF  NM:i:0  MD:Z:47 AS:i:47 XS:i:47

We convert to bam:

samtools view -bS test.sam > test.bam

And read using R:

library(GenomicAlignments)
gr = readGAlignments("test.bam")
gr

GAlignments object with 2 alignments and 0 metadata columns:
      seqnames strand       cigar    qwidth     start       end     width
         <Rle>  <Rle> <character> <integer> <integer> <integer> <integer>
  [1]     chrI      +       42M8S        50       769       810        42
  [2]     chrI      +       3S47M        50      3965      4011        47
          njunc
      <integer>
  [1]         0
  [2]         0
  -------
  seqinfo: 1 sequence from an unspecified genome

The start and the end can be retrieved using:

start(gr)
[1]  769 3965
end(gr)
[1]  810 4011


Is there a good practice (or an easy one for a SBOL newbie) to set up a collection of biological parts and devices? I have a collection of biological parts and devices in a particular format that basically stores information of the sequence of the elements. My intention is to move to a more standardized data structure with a particular interest in Synthetic Biology. SBOL 3 looks like a good choice, then the question is if there is a way to kind of "import" data from my sequence-based database into SBOL?

For SBOL 2, you can readily import data via the online converter tool at https://converter.sbolstandard.org/
The converter basically just wraps the conversion function found in the SBOL libraries, so if you're coding, you can do the conversion directly from the libraries.
For SBOL 3, there isn't yet an online tool, but James McLaughlin has a converter in his JavaScript library: https://github.com/sboltools/sbolgraph
I believe he is working to bring it online shortly, so hopefully this answer will soon need to be updated or replaced by an answer from him.
An anecdote on the difficulty level: I recently started using an SBOL3-->SBOL2-->GenBank conversion using pySBOL3 and pySBOL3. It took a couple of hours to setup the pipeline, and isn't perfect, but it's pretty doable (and ought to be generalized for others to use as well at some point). The other way would be easier, since GenBank contents are pretty simple.


Is the quality score of fastq used somewhere besides trimming/fastqc? In the fastq format every 4k+4-th line contains the positionwise qualityscore (ascii encoded):

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65
...

which program uses this information besides trimming tools (cutadapt) or quality control (fastqc)?

I heard that segemehl (mapping) abandons this quality string in favor of running time (maybe other mapping tools too?)
But what about abyss, trinity, ...?

UPDATE:

So i atrifically generated a fastq file containing 2 identical reads that perfectly match a sequence in some genome (C.fna).
The two reads have the most extreme phred scores (only '!!'=highest sanger phread and only 'II'=lowest sanger phred)

Now using segemehl the mapping quality is identical (60)

@HD VN:1.0
@SQ SN:C_fna    LN:15440
@RG ID:A1   SM:sample1  LB:library1 PU:unit1    PL:illumina
@PG ID:segemehl VN:0.3.4    CL:segemehl -i C2.fna.idx -d C2.fna -q test.fq
A   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
B   0   C_fna   261 60  220=    *   0   0   ACAGGAAACACAGAAAAAAGCCCGCACCTGACAGTGCGGGCTTTTTTTTTCGACCAAAGGTAACGAGGTAACAACCATGCGAGTGTTGAAGTTCGGCGGTACATCAGTGGCAAATGCAGAACGTTTTCTGCGTGTTGCCGATATTCTGGAAAGCAATGCCAGGCAGGGGCAGGTGGCCACCGTCCTCTCTGCCCCCGCCAAAATCACCAACCACCTGGTG    IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII    HI:i:0  NH:i:1  NM:i:0  MD:Z:220    RG:Z:A1 YZ:Z:0
[SEGEMEHL] Thu Feb  6 08:35:15 2020: matching w/ suffixarray has taken 0.000000 seconds.

So segemehl does NOT use the read quality phred (basically only using the read sequence for mapping)

I think almost all alignment programs would use base quality as a factor to calculate mapping quality (See here). These would also be used in Variant Calling to pick reads where high quality bases confidently align as either the reference base or one of the alternate bases.


Creating frequency plot from weblogo with frequency unit on the Y axis Hello I am trying to create frequency plot for a motif with frequency unit on the Y axis from 0 -100 nt as the unit as represented in this paper.(Seiler et al., 2018)

Within weblogo, add the --units probability argument and edit the figure with Illustrator (PDF/SVG) or Inkscape (SVG) to modify the y-axis, changing fractional probability ([0,1]) to a percentage ([0,100]) by multiplying by 100. Refer to the manual for more detail.


Using pythoncyc to get ec-code of a reaction I am trying to collect the ec-code (and gene) of a reaction. I have tried the following:

import pythoncyc
pgdb = pythoncyc.PGDB("sco")
r = "RXN1A0-6128"
gene = "SCP1.233B"
enzyme = pgdb.enzymes_of_reaction(r)
enzyme = pgdb.enzymes_of_gene(gene)
pgdb.full_enzyme_name(enzyme)

However, I only get an enzyme id, and not the ec-code, which is visible on biocyc.

BioGeek is correct, however I will post a full answer.

import pythoncyc
# To get ec-code
db = pythoncyc.select_organism('sco')
reaction = db["RXN1A0-6128"]
ec_code = reaction["ec_number"]

# To get genes
pgdb = pythoncyc.PGDB("sco")
genes = pgdb.genes_of_reaction("RXN1A0-6128")


How many single copy gene are needed to construct the phylogenetic tree? I ran the orthomcl program with genomes of the whole streptomyces genus(from ncbi) as input, but I only got 22 orthologous groups(streptomyces's genome size vary from 2M to 15M, which shows the great difference among them). And the analyst from the bioinformatics company told me that my result wasn't reasonable if I try to construct the phylogenetic tree with the 22 orthologous groups I got. He also told me, I'd better use almost 1000 orthologous groups to build the tree. My question is, how to know my orthologous groups are enough to construct a tree?

The one thing that strikes me is that OrthoMCL was written under David Roos supervision, which really means it was written with Toxoplasma gondii in mind, i.e. a eukaryotic parasite with notable levels of gene duplication. They tested it against Plasmodium falciparum, but that isn't surprising because both are members of the Apicomplexa. Personally I would suggest looking at a tool more orientated towards bacteria, albeit I'd like to see the output of OrthoMCL.

Prokka is a tool that might give better results, here, which is designed for bacteria, but personally I would extract the orthologes/homologues directly if the primary purpose is to contruct a phylogeny, i.e. I would parse annotated genes which are not "hypothetical".


From ensembl id to protein structure Forgive me if this question is too trivial.

I have a list of ~700 Ensembl ID. I need to extract the protein structures of each Ensembl id. How can I do this?

Is there any script(python, R) so that we can download the structure from the website(like RCSB)?

Ensembl to PDB

So the dataset you need is from SIFT. There is one with Ensembl IDs. Majorly this also includes what chain is your peptide and what is the offset. About the latter, I actually wrote a blog post about it recently.

PDB download

Using Python, the Biopython PDB module is very good. Although the PyMol module is more powerful. But as you want only coordinates (right?), you'd best use the former, which has a better documentation.


How to get matching pattern along with ID in a single command in grep? I have a file containing multiple fast sequneces. For a specific consensus pattern as input, I extracted all the matching patterns from target fasta sequences with
grep -o -E "CC[GT]AAA[GC][AC]TT[GC]" input.fasta

However, the above command will retrieve just the matching sequences from fasta sequences and I  also wanted to get the corresponding fasta header of each match.
For example, if input.fasta  file is something like this,
>Gene 1
TGATGAAAAATGATAGAT
ATTGGGGGAAAAAAAAAT

>Gene 2
TTTCCTAAAGATTGT
AAATTTAAAAATGTTTTT

(Gene 2 has matching subsequence CCTAAAGATTG)
Output:
CCTAAAGATTG   Gene2

I prefer a solution with grep. But other possible solutions also helpful.

Use this Perl one-liner:
perl -lne '$id = $1 if /^>(.+)/; ($m) = /(CC[GT]AAA[GC][AC]TT[GC])/; print join "\t", $id, $m if $m;' input.fasta

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-n : Loop over the input one line at a time, assigning it to $_ by default.
-l : Strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches


Genome QC + Assembly Pipeline semantics I’m trying to create a pipeline for genome assembly. How best can I “redirect/pipe” from existing fasta files (or files in general) to other steps of the pipeline?
I was thinking of going from the SRA Download step (prefetch and fasterq-dump) to fastqc, then to trimmomatic or bbduk to remove the adapters, one more fastqc check and then to SPAdes for the assembly.
In my (amateur) mind it would be cool to be like:
“SRRXXXXXX” | prefetch | fasterq-dump | fastqc | bbduk.sh | fastqc (second) | spades.py
Any and all help is greatly appreciate! I know this stuff is pretty basic but I am so new to everything it’s painfully obvious to even me.

Maybe take a look at workflow management systems, like

snakemake
Nextflow
Toil
Cromwell

Admitted, each of those adds their additional learning curve. However, especially snakemake is not so much more complicated than a standard bash workflow. Those systems do help a lot organize and orchestrate your workflows


Manually define clusters in Seurat and determine marker genes I want to define two clusters of cells in my dataset and find marker genes that are specific to one and the other. Is there a way to do this in Seurat? Say, if I produce two subsets by the SubsetData function, is there a way to feed them into some other function that would calculate marker genes? If not, what other packages would you recommend for doing that?

If you look here:

https://satijalab.org/seurat/seurat_clustering_tutorial_part2.html

I just need a way to define ident myself, the number of levels (2) and assign numbers to each cell (0, 1), and then run DE between 0 and 1 clusters which is obvious how to do afterwards.

I think you are looking to FindAllMarkers function from Seurat. As you said, you just have to define your ident, that have to have the structure of a table (cell names as names and cluster as value):

pident=as.factor(clusters)
names(pident)=cellNames
object1@ident=pident

And then run the FindAllMarkers function:

FindAllMarkers(object1, min.pct = 0.25, min.diff.pct = 0.25)

You can specify several parameters in this function (type of DE to perform, thresholds of expression, etc).


HDF5 and BioSQL solutions I'm looking at better database/storage solutions for NCBI virus data, with all attributes particularly year and country of isolation, together with structural data, possible antibody data, T-cell data and bioinformatics data such as %GC etc... The amount of data is comparatively small, with a maximum size of 3 x ([6000 < no. viruses < 20000] for 10 genes) and each entry will contain nucleotide, amino acid and one or more pieces of alignment data (i.e. with indels).

The data will be used for numerous different calculations from structural biology, phylogenetics, and deep learning. Some of the calculations will require parallel processing. A typical calculation would be to identify a clade/similarity/difference then screen the database to look at the isolates information.

Previously I've used lots of flat files which are loaded at initiation and are held within a data pipeline either as an object (Python) or especially a multi-dimensional hash (Perl, preferred).

The two solutions that I think look suitable are:

BioSQL - The official Genbank solution, however my SQL commands are very basic and the documentation appears minimal (for BioSQL)
HDF5 (hierarchical data format) - this looks good because it uses complex hashes (dictionaries), and allows layering of multiple information onto a single Genbank sequence code (key). It has strong documentation. However, I understand this isn't a "database" (long term storage), but thats how I intend to use it.

JSON/Pickling I had a look at JSON, which to me looks okay, because I'm fond of complex hashes / dictionaries, but it appears to be more of a Web solution. Finally, "pickling" (Python) is cool, very cool, but I understand this not a database solution, because it risk compatibility issues, e.g. with future Python releases.

I am interested in knowing what archival systems board members are using, and a critique whether my assessment is ok for moving forward.

Just to conclude, @user172818 and @AlexReynolds are collectively saying either learn SQL or stick with CSV.

Update,
I can provide an answer. It's HDF5 with pandas dataframes (Python).

Simply

df.to_hdf('path/file.h5', 'key', format='table')

This isn't SQL but it works:

What it gives,

it's fast, that's it's purpose (not my main reason though)
it perfectly preserves all pandas dataframe information, this is  important
The 'key' ( object) allows you to store multiple dataframes into the same file, you just provide different key names

The last point is very useful, because you can pull out a dataframe of choice and the SQL commands are mimicked within pandas

I ran into trouble using CSV to store/recover complex pandas dataframe formats.

So it's good enough and much faster than SQL. I could have problems with format changes in time (unlike SQL), but for the here and now it works well.

Just to add a caveat, for dictionaries the json dump works well, but it is just shunting dictionaries around, so isn't relevant.


Histogram of different classes of variables in one column I have a basic ggplot2 R question. I have a dataframe in which data looks like this:

AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7

The file is of multiple lines. I want to make a histogram using gplot2 that will plot the 9 kinds of variables AA,BB,AB,BA,AN,BN,NN,NB,NA on the x-axis and the sum of it's corresponding values in the next column on the y-axis.

So for example, here the y-values should be

AA= 13
BB= 12
AB= 5
BA= 4
NN= 7
AN= 10
NA= 3
BN= 12
NB= 1

This isn't a bioinformatics question, but it's quicker to answer than to close.

df1 <- read.table(text = "
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7", na.strings = "")

library(ggplot2)

ggplot(df1, aes(x = V1, y = V2)) +
  geom_bar(stat = "identity")


How to interpret duplication from MultiQC/FastQC? FastQC of my sample files, aggregated into a single plot by MultiQC. Blue represents unique reads. Black represents duplicate reads. The x-axis is the number of reads.

I see I have more duplicate reads than unique reads, how do I interpret this plot? This is my attempt:

  FastQC took the first 100,000 sequences in each sample (details).
  Generally, higher duplicated read percentage is welcome as it
  represents better sequencing depth and lower sequencing error
  contamination. This dsitribution in this plot looks reasonable and
  nothing unusual.

I would agree with your statement.  Typically the ratio of unique to duplicated reads is a function of input library complexity (number of input molecules) verse sequencing depth (the amount of molecules sequenced).

In the example you gave you have about a 2:1 ratio of duplicates to unique reads which suggests that these libraries have been sequenced to their fullest and additional sequencing would yield few new unique reads.


Transferring genomic features on new coordinates I have a eukaryotic genome for which an updated sequence for a chromosome  was recently obtained. I want to map RNAseq reads on the genome (and perform other downstream analyses) and would like to use the most up-to-date information possible (so the 'new' sequence of that chromosome).

However, before performing the read mapping, I'd like to update the annotation of my genome (GFF3) to match with the 'new' coordinates of this chromosome.

Basically, I would like to transfer information from the 'old-coordinates' GFF3 to obtain a 'new-coordinates' GFF3, and if possible conserve all the information/hierarchy of the file (gene, mRNA, exon, etc.).

Example:

chrXX   source  gene    222 5942    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    222 5942    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    222 5794    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    5889    5942    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 222 5794    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 5889    5942    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

... should be updated to ...

chrXX   source  gene    333 6053    .   -   .   ID=gene_1;Name=gene_1;length=5720
chrXX   source  mRNA    333 6053    .   -   .   ID=gene_1.1;Parent=gene_1;Name=gene_1.1;length=5720
chrXX   source  exon    333 5905    .   -   .   ID=gene_1.1.2;Parent=gene_1.1
chrXX   source  exon    6000    6053    .   -   .   ID=gene_1.1.1;Parent=gene_1.1
chrXX   source  CDS 333 5905    .   -   1   ID=CDS:gene_1.1.2;Parent=gene_1.1;Name=gene_1.1
chrXX   source  CDS 6000    6053    .   -   0   ID=CDS:gene_1.1.1;Parent=gene_1.1;Name=gene_1.1

The approach I tried (but I am not sure it is the way to do it):

Extract sequence of all the features that are going to be updated to fasta (bedtools getfasta)
Map these sequences to the genome with the updated chromosome (gmap with --nosplicing since sequences we are mapping correspond to genomic regions).
Create the updated GFF3 file. When doing that, I apply one supplementary rule: if a feature mapped on another chromosome, but there was also another alignment on the updated chromosome, prioritize the one that correspond to the same (updated) chromosome.

What would be the proper way to do such a thing?  Any suggestion of methods/tools is welcome!

I think the standard way of doing this is to make a chain file, then use it to do a liftOver of the annotations:

mkdir psl
for i in ../ci3/rm/masked/*.masked; do blat ../ci2.2bit $i -tileSize=12 -fastMap -minIdentity=98 psl/`basename $i .fa.masked`.psl -noHead -minScore=100; done

Translate psl files to chains in the directory chain:

mkdir chain
for i in psl/*.psl; do axtChain -linearGap=medium -psl $i ../ci2.2bit ../ci3/ci3.2bit chain/`basename $i .psl`.chain; done

Merge short chains into longer ones into the directory chainMerge:

mkdir chainMerge
chainMergeSort chain/*.chain | chainSplit chainMerge stdin -lump=50

concat and sort the chains:

cat chainMerge/*.chain > all.chain
chainSort all.chain all.sorted.chain

Need info about chromosome sizes for netting:

twoBitInfo ../ci3/ci3.2bit ci3.chromInfo
twoBitInfo ../ci2.2bit ci2.chromInfo

Netting: identify alignable regions from chains:

mkdir net
chainNet all.sorted.chain ci2.chromInfo ci3.chromInfo net/all.net /dev/null

Finally, select the right alignable regions using the nets, creating a "liftOver" file:

netChainSubset net/all.net all.chain ci2ToCi3.liftOver

Run liftOver:

CrossMapy.py bed ci2ToCi3.liftOver test.hg18.bed


How do I efficiently perform a metagenome screen of “all” species? I’ve got an RNA-seq dataset with a large proportion of environmental RNA “contamination”. BLASTing random reads reveals that much of the data comes from bacterial, plant and viral RNA. My target organism only accounts for ~5% of the RNA-seq read data.

I would like to get a more or less (more, if possible) comprehensive overview of what species are found in my sample — bacteria, plants, animals (?), but also viruses. How can I perform this?

I have been unable to find a “standard” way of performing a metagenome screen. BLAST online services all seem to be severely rate limited (certainly unable to handle uploading ~80 M reads). Installing BLAST (or a similar tool) locally is of course not an issue, but I can’t find a comprehensive database that spans all phyla — the best I’ve been able to find are databases, such as NCBI-NR, that are restricted to single phyla or classes, or to bacteria.

Ideally I’d like a workflow that I feed with my RNA-seq data, and get out a list of species with % of reads mapped.

As you point out, BLAST is not going to work well for reads from a typical RNASeq experiment.

From an RNASeq dataset, you might be able to get away with just a ribosomal screen using Silva, because the ribosome is the most abundantly transcribed gene. This will usually be the case even if polyA selection or ribosomal depletion was carried out on a sample.

However, bear in mind that there will be a lot of false results due to database incompleteness. I used Kraken when I did something like this a few years ago, and ended up with one sequence annotated as a fruit bat (where my subject was a water-living planarian flatworm). Centrifuge would be a better option now; the authors even claim that they can store the entirety of nr in a 70 GB index with Centrifuge (which would be better / more comprehensive than just a ribosomal screen).


How to obtain clusters of hierarchical heatmap when using Python? Is there a good way of obtaining the labels (e.g. genes) within individual clusters that haven been clustered hierarchically in Python (preferentially, but not necessarily, by seaborn)?

I found these difficulties highlighted elsewhere (here, and here), and would hope for some simple option (e.g.: like in MATLAB where cutoff can be specified, and corresponding genes would be returned), without potential risk that the plotted hierarchical heatmap differed from the cluster assignment ( see here)

Clustering like this is typically done with scipy. Here's the code we use in deepTools (original context here):

from scipy.cluster.hierarchy import fcluster, linkage
Z = linkage(some_matrix)  # You might want to set `method` and `metric`
groups = fcluster(Z, nGroups)  # You might want to set `criterion`

groups is then a vector containing the group assignment of the input matrix according to the clustering. As an aside, this is exactly the same as what's shown in the last link in your post, so I'm not sure how you didn't find it.

BTW, at least in python, hierarchical clustering is painfully slow with large datasets.


Is it possible to go from a digital sequence to an actual nucleotide? I'm trying to learn about what the state of technology is at presently.   It seems like we clearly can go from nucleotide to digitally stored sequence, but can we transcribe something from the digital sequence?

Oligonucleotide synthesis is a service offered by many providers such as IDT.


using SNPs to identify mixed samples Is there a way to identify mixed samples based on SNPs?

Example input (table of genotypes for multiple samples):

|      | sample1 | sample2 | sample3 |
|------|---------|---------|---------|
| rs1  | AA      | BB      | AB      |
| rs2  | AA      | BB      | AB      |
| rs3  | AA      | BB      | AB      |
| rs4  | AA      | BB      | AB      |

Desired output: Sample relationships/distances. In this example, sample1 and sample2 are distinct, but sample3 is likely a mix of two (biological or technical).

Based on this simple example, the calculation seems doable, but maybe it's impossible in more complex realistic scenarios.

There is a nice tool peddy that does something conceptually similar, but it expects pedigree information as well as genotype frequencies for both the input samples and population.

Yes. The general term for identifying subsets from a mixed group is deconvolution.

This is easier to do when there is haplotype-level information available (e.g. sequenced reads that span multiple SNPs), but it is possible to get some indication of sample mixes from unlinked SNPs alone.

Deconvolution is helped a lot by having proportional data in the input (e.g. frequencies of genotypes, rather than a single genotype). Most SNPchip platforms have multiple probes per variant, and can provide some indication of variant proportion when digging down into the raw data. Deconvolution is also a lot easier with some knowledge about the component mixes. It's much easier to simulate mixes than to guess at possible mixes.

Here's a paper about the application of deconvolution in forensics.

Here's a paper about a program called CIBERSORT which carries out a deconvolution of immune populations based on gene counts:

Here's another paper which discusses the idea in more general terms.

You might find more information on one of these Biostars posts:

Deconvolution Methods on RNA-Seq Data
Differential gene expession analysis in cell populations of mixed tumor and normal cells


How to efficiently get human gene names from NCBI based on a large list of SNPs I found a good answer related to my question here:
How to get a list of genes corresponding to the list of SNPs (rs ids)?

But it says about small number of SNPs. I want to get gene names based on approximately 12000 to 19000 SNPs. NCBI says about 10 read per second with API key. Is there any efficient way to get gene names for thousands of SNPs? Is that database available to store locally? Any advice please?

Modify SNPs and gene annotations for your genome assembly of interest. This example shows how one might do this for hg19, for instance:

Download SNPs and convert to BED:

$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz \
    | gunzip -c \
    | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - \
    > hg19.snp151.bed

Filter this BED file for entries of interest:

$ grep -wFf snps-of-interest.txt hg19.snp151.bed | cut -f1-6 > hg19.snp151.filtered.bed6

Download gene annotations and convert to BED:

$ wget -O - ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz \
    | gunzip -c \
    | awk '($3 == "gene")' \
    | gtf2bed \
    > gencode.v19.genes.bed

Note: In the future, newer versions of Gencode annotations may be available for your reference assembly. Always check their site to confirm what is most recent.
Map genes to SNPs:

$ bedmap --echo --echo-map-id hg19.snp151.filtered.bed6 gencode.v19.genes.bed > snps_with_associated_gene_names.bed

A list of genes that overlap a SNP of interest will be in the seventh column of the output.

If this must be done within Python, you could use something like subprocess.call(...), where ... includes commands wrapped in bash -c.

I don't recommend this because of the need to escape quote marks and deal with Python's quirks in running command-line tasks, which makes a Python-based approach very difficult to set up and maintain. I'd suggest learning some basic bash scripting to make this easy and fast.


What is a PacBio "movie file"? I came across references to "movie files" from PacBio sequencing in this paper:
https://www.jimmunol.org/content/204/12/3434
Specifically:

Movie files used to generate results presented in this article have been submitted to the National Center for Biotechnology Information BioProject (https://www.ncbi.nlm.nih.gov/bioproject/) under accession number PRJNA389440.

Looking at the SRA entries, it looks like it's just the raw reads as I expected.  But what is meant by "movie files"? (This is probably just my ignorance about PacBio as opposed to Illumina or others but I've never seen that term before and can't find much online.)

[disclaimer: not a PacBio person, I'm just vaguely aware of the technology]
Movie files are the PacBio raw data format, representing observed signal intensities in each sequencing well over the course of the sequencing run. The movie files can be used to re-call sequences with updated statistical models, which can improve accuracy, or test out hypotheses about how epigenetics alters the rate of nucleotide incorporation (i.e. methylation).
Illumina has a similar static version of this which is an image (or multiple images, depending on chemistry) of the flow cell after each synthesis cycle. Illumina's raw data format consumes tens of terabytes per run, and is of minimal value after base calling, so is almost always discarded after base calling.
Anything uploaded to NCBI SRA must contain FASTQ information. While the original uploaded files can be kept (which may include time/signal information), the files are converted into a FASTQ-like format for NCBI's internal representation. As such, there's a chance that movie files could be uploaded to SRA, but those files would also need to include the called FASTQ information (e.g. in a HDF5 container). In contrast to this, ENA does allow uploading raw data files without called information.


Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python) I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.

I have a population allele count data, and it looks like this:

1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0

Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").

In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.

I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:

1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0

As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?

Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...

awk '{for(i = 1; i <= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i<=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i <= NF; i++)
        allele2[i] = $i;
            for(i = 1; i <= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i <= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'

I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:

allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000),
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

Edit

I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution

I used readr for file IO:

library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

Here are the timings:

cat Benchmarks/R.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.


convert tRNA coordinates Result to Fasta format I believe many of the experts here have gone through this issue,

I have run tRNAscan-SE on my genome to predict all the tRNAs, and have successfully gotten all the tRNA coordinates like start, end and type all  arranged in a table.

My question is how can I convert this result to fasta format as long I have all the details about the tRNA ready?

Is there any tool that I can use to get the fasta sequences of the regions described in my table from the original genome?

Sequence                tRNA    Bounds  tRNA    Anti    Intron Bounds   Inf
Name            tRNA #  Begin   End     Type    Codon   Begin   End     Score   Note
--------        ------  -----   ------  ----    -----   -----   ----    ------  ------
NW_011590950.1  1       36704   36632   Ala     AGC     0       0       40.8    pseudo

You can select “Output BED format” on the tRNAscan-SE server:

Given the output BED format and your input reference in FASTA format, you can easily generate a FASTA file of tRNA sequences via bedtools getfasta:

bedtools getfasta -name -s -fo ‹output.fasta› -fi ‹reference.fasta› -bed ‹tRNAscan-SE-output.bed›


How to find the sum of all sequences in a multifasta file I have a multifasta file containing several genes from a genomic island. I want to know the total length of this island. How can I do this?

grep -v ">" multi.fasta | tr -d '\n' | wc -c

grep -v ">" find all non header lines (i.e. without ">")
tr -d '\n' delete all newlines
wc -c count all characters


Common gene naming conventions and how to convert between them I have two questions, one is a direct question about some RNA expression data.  The other is more broad, but motivated by the first.  I have downloaded an expression set for a class project and I am unfamiliar with the naming convention, I recognize entrez ID but not this one.
    > rownames(Rao)[1:10]
 [1] "0610007P14Rik" "0610009B22Rik" "0610009D07Rik" "0610009L18Rik" "0610009O20Rik" "0610010B08Rik" "0610010F05Rik"
 [8] "0610010K14Rik" "0610011F06Rik" "0610012G03Rik"

So my first question is what is this gene naming convention.  The second is what are many other conventions I should recognize and what tools can I use for converting between them?  Lastly what kind of errors may crop up when I do that and how can I check for mismappings or non-unique mappings?
I am primarily interested in R but methods in python would also be appreciated!  Its always good to have more tools.

Biologists are familiar with SYMBOL (for good reasons).
Bioinformaticians like ENTREZID ENSEMBL REFSEQ GENBANK etc, I'm definitely missing some here.
You don't need to know these by heart, just know that databases don't map 1:1.
I've worked with mouse RNAseq data before, so I agree with user::RamRS comment that it looks like a mouse gene. What type of identifier is it? I had no idea.
You can find this out at NCBI::GENE, by querying for, let's say Rao[1] "0610007P14Rik". (https://www.ncbi.nlm.nih.gov/gene/?term=0610007P14Rik)
NCBI::GENE tells you the SYMBOL::Erg28, ENTREZID::8520, and in the far right a column called ALIAS::0610007P14Rik.
Now that you know your identifier type (ALIAS), you can map between databases in R.
There are many packages to do this. Someone mentioned biomart in the comment.
I use AnnotationDbi
library(AnnotationDbi)
library(org.Mm.eg.db)

Rao <- c("0610007P14Rik", "0610009B22Rik", "0610009D07Rik", "0610009L18Rik", "0610009O20Rik", "0610010B08Rik", "0610010F05Rik", "0610010K14Rik", "0610011F06Rik", "0610012G03Rik")

#use the function select to retrieve the column type "SYMBOL" with keys "RAO".
#"RAO" is a keytype of "ALIAS". And I want to search in the mice database "org.Mm.eg.db"

AnnotationDbi::select(x = org.Mm.eg.db,
+                       keys = Rao,
+                       keytype = "ALIAS",
+                       columns = "SYMBOL")
'select()' returned 1:1 mapping between keys and columns
           ALIAS        SYMBOL
1  0610007P14Rik         Erg28
2  0610009B22Rik 0610009B22Rik
3  0610009D07Rik         Sf3b6
4  0610009L18Rik 0610009L18Rik
5  0610009O20Rik         Dele1
6  0610010B08Rik        Gm2026
7  0610010F05Rik 0610010F05Rik
8  0610010K14Rik 0610010K14Rik
9  0610011F06Rik       Mettl26
10 0610012G03Rik 0610012G03Rik

In the cases when they don't map 1:1, the function will throw a warning, and you can access it by warnings(). It is rarely important though.
Hope this answers the question.


Error in createMAE function: non-unique values when setting 'row.names' in TCGA LIHC Data I am working with the createMAE() function of the ELMER Bioconductor package. While executing the createMAE() function using TCGA LIHC DNA Methylation 450K and RNASeq HTSeq - Counts data, I am getting this error:

Error in.rowNamesDF <-(x, value = value) :

duplicate 'row.names' are not allowed

In addition: Warning message:

non-unique values when setting 'row.names': 'TCGA-DD-AACA-02A', 'TCGA-DD-AACA-02B'

Can anyone tell me how to get rid of such error?

The error message is pretty clear: R requires unique row names on data frames and some of the "values" that you are trying to set as row names are non-unique.

You will need to provide a vector i) of length the same as your row number, ii) that is composed of non-unique values.


IGV faster loading remote bigWig files I noticed remote bigWig files (http URL) take sometimes about 1 minute to load on IGV. Once they loaded, jumping from one region to another can also take several seconds.

My understanding is that bigWig files have an internal index for different chunks of different sizes. Can the index be modified so that the speed of loading on IGV is improved? Does this need to happen at the point of creation of the bigWig file, or can the index of an existing bigWig be updated to speed up loading?

You would need to do this when creating the bigWig file. Note that I'm not aware of anything that allows manually setting the various zoom levels (libBigWig/pyBigWig will allow you to set a max number, but that's it).

This then raises the question of whether IGV is using an appropriate zoom level to begin with. I've never looked at how it chooses which of the zoom levels to use, but I know it uses them since if you don't have any it will fail to load bigWig files.

Honestly, though, if it's taking a minute to jump to a new region then I suspect that the network is just slow or the remote server is overtaxed. Changing the internal indexing isn't going to fix that.


Survival analysis using CoxPH - Effect of covariates Hi and sorry for the long post in advance,
I'm doing a survival analysis of lung cancer patients using Python's lifelines package.
According to the documentation, the function plot_partial_effects_on_outcome() plots the effect of a covariate on the observer's survival. I do understand that the CoxPH model assumes that the log-hazard of an individual is modelled by a  linear function of their covariates, however, in some cases the effect of these covariates seems to perfect when plotted. I'll illustrate these with three images below describing the issues I have in my data analysis.
Example 1:
This one I'm providing as it looks realistic to me, there is a clear difference in the effect of the covariate when positive or negative.

Example 2:
This one just looks gimmicky, my concern is firstly there seems to be a perfect relation between the pathologic stage of cancer and the survival, it's logical that higher stages should impact the survival negatively, however, I don't understand why the graph is so symmetric as the survivability of cancer shouldn't be proportional to the patient cancer staging.

Example 3:
Same problems as for example 2, this time with a binary not categorical variable (so I'm less reluctant on disbelieving it), however, the curves do seem to symmetric, almost like it's the same curve but with a standard deviation. I also highlighted a region of interest which I don't understand why both curves stabilize (maybe lack of patients in those regions?)

Example 2 and 3 seem symmetric because they are forced to by the model. The Cox proportional hazard relies on the proportional hazards assumption. This basically means that the ratio of the hazards are constant over time. If you were to plot the hazard ratio, it would look like a flat line over all Event_at. Another way of saying this is that the hazard ratios are constant over time.
How that shows up in a survival plot is a little different because of the relationship between hazards and survival. The survival curves end up following the pattern exemplified in Example 2.
As for example 3, I would guess that there are few individuals with event times past 3500. This is why you have those large jumps.
Proportional Hazards
So what does the proportional hazards assumption mean? Below are some plots of hazard functions to help your intuition on what the proportional hazards assumption means. In all of the following scenarios, the hazard ratio is 2 for blue-to-red. However, the proportional hazards assumption is only true in A, B, and C.

Since the proportion is remaining constant, it can look a little hard to see. To make it easier on the eyes, you can log-transform the hazard functions. Below is a plot that does that.

As you can see, the hazards are a constant distance from each other for A, B, and C. Basically, the proportional hazards assumption is saying that your scenario looks something like one of those plots.
In scenario D, the hazard ratio varies over time. At the beginning it is basically HR=1. At the end it is basically HR=4. With the proportional hazards assumption, Cox regression will essentially take the 'average' of the hazards over time. In this case it would give you an estimate of HR=2. Depending on the scenario, this 'smoothing' may be more or less acceptable. But this is basically what the proportional hazards assumption does in the background.
The proportional hazards assumption is actually more problematic in the cases where the two hazard functions cross each other. Below is an example of that.

In this scenario, blue has an increased hazard at the start (relative to red) then a decreased at the end. A Cox proportional hazards model would give you HR=1 in this case, because the average over all time points is no relationship.
So how do you tell if you violate the proportional hazards assumption? You may want to use a stratified Kaplan-Meier by your variables to get an idea of what the hazards look like without stipulating the proportional hazards assumption. You can use the Kaplan-Meier estimator to get a non-parametric estimate of the discrete hazard (i.e. it does not make the proportional hazards assumption). Then you can log-transform those hazards to get a plot like the second. There are also statistical tests and other diagnostics (Schoenfeld residuals) to check this assumption.
The next question is do those violations matter? This somewhat depends on the question you are trying to answer and how substantial the violations are. In cases where the two log-hazard curves don't cross (like in D), it may be acceptable to smooth over the variations. However, it settings where the curves cross in the final figure, it is probably best to not ignore this problem.
The original write-up of this answer was here


Evaluate clusters of individuals by using their sequence data For a dataset of several hundred individuals, I applied a hierarchical clustering to generate clusters based on a functional trait that sets them apart.

My task is now to evaluate if these clusters can be supported by the nucleotide  sequence data of the corresponding gene, i.e. is there more genetic similarity within each cluster than between the clusters.

For a first approach, I created a multiple sequence alignment for every cluster and calculated the the % identical sites value. There are more identical sites within each cluster than in an alignment of all sequences.

Do you know any good tools (preferably python) that can perform a more sophisticated evaluation of existing clusters based on sequence data?
Are there other scoring methods with which I can evaluate the sequence similarity within the clusters (again preferably python)?

You can calculate allele frequencies for each cluster you have to further verify if they belong to similar population, however if size of your dataset is rather small this may not work for you.

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3813878/, here is an article about differentiating populations based on sequencing data.


Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR? I'm writing a python function from scratch to do this, but I feel like this must exist in some standard bioinformatics library already. In principle, this is a simply regex operation which many must have written previously. With the goal of having a centralized reference as opposed to everyone writing their own scripts to accomplish the same task (which is time-consuming and not necessarily reproducible), here is the "construct a CIGAR" function:

Here's an alignment:

ACGT-TGC

Here is a reference:

ACGTATGC

The reference will be required in order to record indels. One could in principle have a series of "queries", e.g. an R vector or a Python list.

This would generate the following CIGAR:

4=1D3=

Which R packages/ python library exists to do this operation? I suspect there is something related to Bioconductor as well.

Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:

def gen_cigar(ref, qry):
    if len(ref) != len(qry):
        raise Exception('unequal length')
    cigar = []
    for i in range(len(ref)):
        r, q = ref[i], qry[i];
        if r is '-' and q is '-':
            raise Exception('both gaps')
        op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
        if len(cigar) > 0 and cigar[-1][1] is op: # add to the last operation
            cigar[-1][0] += 1
        else: cigar.append([1, op]);              # a new operation
    return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string

print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));

The output is:

4=1D1=1I2X

PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.


Is there an established method for comparing expression of groups of genes (gene sets)? I have a figure in a paper where I show the log2(fold-change) (obtained with DESeq2) between two groups (based on genotype) for genes within a specific hallmark gene set, and a reviewer is asking whether the differences are statistically significant.
I have not been able to find any methods for comparing expression of groups of genes like this, and i'm sure there is a method or a reason why it isn't done that I'm not aware of.
To be extra clear, what i have is two sets of log2FC values for a set of genes (one for one genotype (drug-treated vs. untreated) and the same for a second genotype), and I'm looking for a statistical comparison of these two groups of Log2FC measurements.
Thanks

You could apply the GSVA transformation to later compare them. You would need to give a step back and start from the expression values not log2FC, but I think it would be better than using the log2FC directly.
If you need to use them you could compare the gene set enrichment value of these two sets on the same comparison. I would use fgsea with the log2FC and the two sets of genes. If the normalized enrichment value (NES) is different then you could explain that these sets expression are indeed different.


Generating random BED intervals given constraints Problem is to generate a random BED interval given the following constraints:

minimum start
maximum end
fixed length
maximum number of masked bases (similar to -maxN option in faSplit)
set of intervals to avoid overlap with
stay within chromosome size constraints

I can think of computationally intensive ways to code it up from scratch, but am wondering if there is a more efficient approach before re-inventing the wheel.

seq_records = {x.name: x for x in SeqIO.parse('path/to/genome.fa', 'fasta')}

def generate_random_interval(chrom, lower, upper, length, maxrep=1.0, avoid_intervals=None):
    while True:
        start = np.random.randint(lower, upper - length)
        end = start + length
        regenerate = False

        for interval in avoid_intervals:
            if (interval.start < end) or (interval.end > start):
                regenerate = True
        if 1. * seq_records[chrom].seq[lower, upper].count('N') / length > maxrep:
            regenerate = True

        if not regenerate:
            break

    return chrom, start, end

You can perhaps collapse requirements for (1), (2), and (6), so long as bounds from (1) and (2) fall within the bounds of (6).

To help with this, you could use fetchChromSizes from the Kent UCSC utilities to quickly get chromosomes and chromosomal bounds for your genome of interest.

Given a chromosome and a maximumEnd that falls within its bounds, uniformly sample from [0, maximumEnd-minimumStart-maxLength), adding minimumStart to determine start position and maxLength as length to get the stop position.

To obtain the desired number of samples in an efficient way, you could pipe these "candidate elements" into bedmap --echo-map-size operations on repeatmasked regions (4) and bedops -n 1 on "regions-to-avoid" (5) operations to do rejection sampling.

Perhaps something like:

$ REPEATMASK_THRESHOLD=123
$ SAMPLES=1000
$ someScriptThatGeneratesCandidateIntervalsWithinBounds \
    | sort-bed - \
    | bedmap --echo-map-size --echo --delim '\t' - repeatmaskedRegions.bed \
    | awk -vthreshold=${REPEATMASK_THRESHOLD} '($1<threshold)' \
    | cut -f2- \
    | bedops -n 1 - regionsToAvoid.bed \
    | head -${SAMPLES} \
    > qualifyingCandidates.bed

This will be relatively efficient, because once the head process at the end writes the number of SAMPLES to qualifyingCandidates.bed, it will send a SIGINT signal up the pipeline that will terminate the upstream processes. So you're only doing sampling to the extent needed.

Hope this helps!


RNA_Seq Analysis in R, propmapped function issue I'm currently trying to learn RNA-seq data analysis and differential expression in R. I've been using the Combine lab tutorial.
Everything was running well until the tutorial asked to use the propmapped function (on this page) on my BAM files and then view the promapped proportion for each sample. In this tutorial we should be getting promapped proportion ranging from 0.05 to 0.150, however in my case I have no reads that align to my mapped genome (promapped = 0 in my case, for each sample). I don't know why as I've been using the data they provided for this tutorial.
I was wondering if someone could follow the first few steps of this tutorial until the propmapped step and tell me if they eventually get the right proportion of reads or not after the propmapped function.
This could tell me if the tutorial is the culprit or not.

There must be a problem in the BAM files you are using. Here are steps I would do to figure where exactly is the problem.

Check the data you have downloaded. The tutorial is using 12 fastq.gz files (SRR1552444 - SRR1552455). You can check on sra how big they are supposed to be (For example SRR1552444 is supposed to be an ~1.7Gb file) and if your downloaded files are about the same. Check if the files are fine (In the terminal you can type zcat data/SRR1552444.fastq.gz | head to check the first few reads). If the files are 0 sized or somehow corrupted, you need to redownload them.

Check if your mapping reference looks alright (size wise, look at the first few lines) and if the reference is indexed (there is a .fai in the directory).

Check the size of the BAM files, they also need to be in the range of a few hundred Mb to a few Gb. I expect these will be corrupted for sure, otherwise, you would get your numbers. You might also want to check if they are indexed (although this should not be a problem, I suspect they will get indexed automatically). If there is nothing wrong the bamfiles, you might want to check another tool to summarise a bamfile (like samtools flagstat) to check if at the end the problem really is in propmapped. If that's the case.

This should at least help you to zoom to the core of your problem if not resolve it.


Why doesn't an (Entrez eutils) einfo request for "gene" return the link gene_nucleotide or gene_nucleotide_pos links? I'm updating a galaxy tool wrapper for Entrez's eutils suite and I'm trying to create a form with valid link selections (among other things) based on the "from" & "to" databases to reduce the size of the select list.
Einfo doesn't return these links for the gene or nuccore databases. I'd found this resource which seemed to show all the possible links, but it it doesn't show these links either. Yet, in working on elink's acheck feature, using some gene IDs in a test run, I saw 2 links that were not listed in either resource:

gene_nucleotide
gene_nucleotide_pos

I manually checked that I can use those as link names in valid elink queries. It did return valid results (i.e. no error).
So I ran elink on the gene database using a few random UIDs to see if I get back a comprehensive list of links (albeit possibly empty results for those links), but it also was missing these 2 links.  I'm concerned that I'm not getting all of the possible links for my form interface and that it could produce valid results and that users may want to link via these links, given prior knowledge that they exist.
How do I get a comprehensive list of database links from Entrez and why doesn't einfo return all possible links?

I've been corresponding with NLM about this issue and I finally took the time to try out their suggestion (which personally I found hard to see between the lines and which is not a discrete solution, but rather a very time-consuming manual process containing false positives because they say that to get a formal and comprehensive response to my query:

Collecting the linkname will be a difficult task, that will take time and coordinate/check with relevant groups maintaining individual databases. Your patience will be greatly apprecired. [sic]

It took me a bit to decipher their suggestion on how to manually find all links, so I will share what I've learned.  There is (currently) no codified or formal means to obtain all possible link names and I infer that many links exist simply because they are utilized in the internal workings of their website.  I could be wrong about that, but suffice it to say that there are many undocumented links and they change constantly.
You can get a list of filter items using each individual database's advanced search web interface, "most of which should be the linkname from that source database to other target database".  (So this should result in all possible links among a series of false positives.)
Here's how you do that.  Let's take the gene database as an example:

Go to the advanced search web page for that database, e.g. https://www.ncbi.nlm.nih.gov/gene/advanced
In the "Builder" section, select "Filter"
Click "Show Index List" to the right of the search term field next to the select list where you selected "Filter" (you will see a multi-select list appear below the search field)
Click "Next 200" until you find a selection in that new select list matching the database name followed by "all" and a number in parentheses, e.g. gene all (29994947) and click to select that item
Click "Refresh index" at the lower right corner next to the multi-select list

The multi-select list will be repopulated and most of the items will be linknames (though you have to replace spaces with underscores in order for those linknames to work with the elink utility).
Doing this for the gene database and then scrolling through just links to nuccore (and its alias "nucleotide" - ignoring est and gss), you will find:

gene nuccore
gene nuccore clust
gene nuccore mgc
gene nuccore pos
gene nuccore refseqgene
gene nuccore refseqrna
gene nucleotide
gene nucleotide clust
gene nucleotide mgc
gene nucleotide mgc url
gene nucleotide pos

If you look only at the linknames in the documentation or via einfo (by supplying the gene database and looking at similar links to nuccore), you only get these links:

gene_nuccore
gene_nuccore_mgc
gene_nuccore_refseqgene
gene_nuccore_refseqrna
gene_nucleotide_mgc_url

So I believe that given this information and the ever-changing nature of these links, I believe I will allow the user the option to enter a linkname manually, if the set in the select list I generate via einfo does not contain the link they need.


how do I convert a vcf file for a GWAS study (using R package vcfR) I have a vcf file with individual level genetic data, that I read using the R package vcfR

gen <- read.vcfR('test.vcf')

Scanning file to determine attributes.
File attributes:
  meta lines: 9
  header_line: 10
  variant count: 152
  column count: 3798
Meta line 9 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 152
  Character matrix gt cols: 3798
  skip: 0
  nrows: 153
  row_num: 0
Processed variant: 153
All variants processed

How do I convert this file to a data frame that I can use for a GWAS or Mendelian Randomization study?

I was told that this should be a data frame where each column is a genetic variant and it take values 0, 1, or 2, depending on the number of "mutations".

I never worked with genetic data before and do not know much about genetics.

If you want to do GWAS or mendelian randomisation, you can do it with Plink (v2) which should be faster than R.


How to get unique somatic mutations for each individual patients I have exome data 32 patients. My aim is to identify unique somatic mutations in each patients. For this I have completed my pipeline for Whole Exome Data analysis using various tools like BWA mem, picard, gatk etc. Now I have VCF and gVCF file of each patient.

For Variant Quality Score recalibration (VQSR), as per the doccumentation, we should have at least 30 patient data for good recalibration. So in order to perform VQSR, I combined gVCF of all 32 patients using genotypeGVCFs tool of gatk and done with VQSR. But now I have cohort of all patient data (a large combined vcf file).

I don't know how variant sub-class tool or other extraction tool work to extract all of the unique and related features to each patient. I know how to get significant mutations using mutect but here it is combined file but I want to get it for each one individually.

How can I get a vcf file for each patient (without missing a single variant associated to that patient) ?

If what you want is to split the main VCF file into 1 file per sample, you could use bcftools query and view commands.

A similar question was asked on biostars, adapting Jorge Amigo's solution to your situation, you might use:

for sample in $(bcftools query -l $combined_vcf)
do
    bcftools view -c1 -Ov -s $sample -o ${$sample.vcf} $file
done


Extracting p-values using corrr package in R I'm performing correlation analysis using the corrr package in R, but I'm unable to extract the p values. Does anyone know how to get the p values in corrr package in R?

I would suggest looking at some of the documentation.
As far as I can tell, the intention of the package is to visualize and organize correlation coefficient estimates. It does not seem that the author of the package has any interest in p-value estimation.
If you are interested in computing p-values, I would suggest using a different tool, such as the builtin cor.test() function.


samtools mpileup empty when filtering out flags I produced a bam file by aligning reads to a small set of synthetic sequences using bwa-mem.

I am heavily filtering reads that are not paired and of a certain orientation.
Applying the filtering, I get a few thousands of reads:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | wc -l

Gives me 89502 reads.

If I then pipe this into samtools mpileup, I get no results:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools view -h -f32 -  | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 0.

I tried different combinations of filtering, and when I do both -f 16 and -f 32 returns empty, but if I do either of those, then it works:

samtools view -h $myfilebam | \
samtools view -h -F4 - | \
samtools view -h -F8 - | \
samtools view -h -F256 - | \
samtools view -h -F512 - | \
samtools view -h -F1024 - | \
samtools view -h -F2048 - | \
samtools view -h -f16 - | \
samtools mpileup --excl-flags 0 -Q0 -B -d 999999 - | wc -l

Returns 1056.

Any ideas why? My thinking was that it would work with --excl-flags 0.

EDIT: substituting mpileup for depth does work, and prints out each position and the depth as expected.

EDIT2: adding -q 0 to mpileup gives the same empty result.

Thanks in advance

By using -h in the samtools view command, you're including all the header lines in your word count. If you happen to have about 89500 reference sequences, then the lengths of those would all appear in the header and inflate the -h word count, but not the mpileup count. Try piping it through an additional samtools view (i.e. without -h) and see if the counts change:

...
samtools view -h -f32 -  | \
samtools view | wc -l

Also, samtools mpileup by default only considers high-quality bases and concordant reads. Try adding a -A to your mpileup line (which stops anomalous read pairs from being discarded):

...
samtools mpileup -A -Q0 -B -d 999999 - | wc -l


How to deal with correlated SNPs in GWAS? I have population SNP data in VCF format, and I found that some SNPs have a great similarity across samples(> 99%). For example:

CHROM_POS   s0  s1  s2  s3  s4  s5  s6  s7  s8  s9
chr1_1  A   G   G   A   G   G   A   A   G   A
chr1_2  A   G   G   A   G   G   A   A   G   T
chr1_3  C   C   C   A   C   A   C   C   C   C

The similarity between chr1_1 and chr1_2 is 0.9, because they have a different SNP in sample s9.

Is there a good way to remove these similar SNPs before I put the data into another pipeline?

Added: These SNPs will be used for GWAS analysis.

What you are attempting to do is known as LD-pruning.
As @Emily_Ensembl said, it is not customary to do this for standard association tests: it is possible that one of the SNPs you remove is causal, or a better proxy for the causal locus, and would give (slightly) better association signal than the other. Even for SNPs in perfect LD, pruning is unwise because it would complicate interpretation: if SNP chr1_2 had some obvious function, but you removed it and are left only with its proxy chr1_1 (which could be far away), you will have a hard time connecting the dots. Typical GWAS tests are handled pretty quickly on modern computers, and are linear in SNPs, so there is no computational reason to prune either.

On the other hand, many fancier procedures would require LD-pruned data, because they either
a) are very computationally expensive, or
b) assume independent SNPs.

I suggest reading PLINK's LD pruning documentation, as it pretty much covers the standard ways to do this.


Removing genes with less than a correlation cut-off between two matrices I have two matrices like this:

> head(DEG_log[,1:4])
            cl1       cl1       cl1       cl1
ALB    4.653796  7.046317  7.732135  7.758708
AQP9   7.798208  8.231112  8.734799  8.054087
CALML5 5.791524  9.058341  6.977369  9.503477
CCL4   9.739961  7.663303  8.934939  7.324958
CDK6   9.832115  9.710043 11.439345 10.330253
CHGA   5.139340 12.024606  5.988773  8.696836
>

I have the same patients and genes in two separate matrices coming from two sequencing panels with the same chemistry. Each sample has a match in another panel. I want to remove genes that have low correlation between two matched samples.

I have normalised and log-transformed the data by DESeq2.

I just plotted the correlation between each pairs of samples:

library("ggpubr")
ggscatter(my_data, x = "sample1", y = "sample1_1",
          add = "reg.line", conf.int = TRUE,
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "onco", ylab = "biomarker)")

Like this:

by:

corr_test(a, b, 0.95)

I can calculate correlation, but I don't know how to remove genes with low or non-significant correlation between matched samples.

Nice data and good question.

It is called "outlier" analysis or the analysis of the residual. Its an easy analysis, but . You calculated the standard deviation using the regression slope as the 'mean' and retain data within 1.96 standard deviations. You can plot the data like this and if you are lucky will see the residual forming a nice 'normal distribution'. Outliers outside 1.96 sd can be formally excluded.

Calling (and plotting) the residual from the regression slope will be possible in R, but I use external packages because its easy. BTW you should also try Spearman's coefficient to check you are getting the same answer.


Extract residue sequence from pdb file (in biopython, but open to recommendations) I'm new to Biopython and I'd like to extract the sequence of residues from a pdb file.
My two questions are:

What is the simplest way to do this? (Esp. when there is more than one sequence)
and
Should I be worried about "PDBConstructionWarning: WARNING: Chain B is discontinuous"?

So far, I've obtained the residue sequence via:
p = PDBParser()
structure = p.get_structure("1ppi", "1ppi.pdb")
ppb=PPBuilder()
for pp in ppb.build_peptides(structure):
     print(pp.get_sequence())

seq = pp.get_sequence().__str__()

This seems to work well for this molecule. However, is there an easier way, especially when there is more than one sequence?
For example, I've read that one can also do
res_list = Bio.PDB.Selection.unfold_entities(structure, 'R')

but res_list is not a sequence of residues in str, and I don't know how to convert the output from res_list to a str sequence.
In addition (or perhaps because I'm going through the PPBuilder), I've recently gotten a lot of warnings of the type:
/usr/local/lib/python3.8/site-packages/Bio/PDB/StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line..
(For example, with structure = p.get_structure("5owu", "5owu.pdb"))
I've seen a discussion about how to silence such warnings, but should I be worried about this? I've also noticed that when I get these warnings, pp builder seems to give me more sequences then are there.

There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

Try help(PDBParser) for more, or see:
http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#__init__


Variant calling for a subset of genes using whole gneome sequencing data I have few 100 raw fastq files from whole-genome sequencing data and I would like to map these files to a set of genes only (and not whole genome) so as to find SNP's associated with them. Can anyone tell me what snp analysis pipeline is best for targeted variation and not whole genome?
I have this in mind:
Mapping the reads with bwa and then using varscan
Hope to hear for more options

Your approach is not recommended. Always align against the entire genome. If you align against a subset you might accumulate false mappings. The aligner will always try to find the best match for every read. If the true origin of the read is not in the reference then the aligner will still try to find an acceptable mapping position. Therefore, align against the full reference and then extract the regions you are interested in, e.g. with the region option of samtools view:

samtools view [options] <in.bam>|<in.sam>|<in.cram> [region ...]

Check its manual for details.

While there is basically nothign wrong with VarScan it is quite old and not maintained anymore. I suggest you give bcftools a try (check its manual). Other maintained alternatives can be (among others) freebayes, strelka2, VarDict or the GATK from the Broad institute.


Interpreting Intergrative Genomic Viewer (IGV) I was following a tutorial on "Tuxedo Genome Guided Transcriptome Assembly Workshop" and was wondering how to interpret the following:

From what I understand from 'Color Legends', the color blue represents something that is below normal and red represents something that is above normal.

Also, my understanding is that merged.gtf is the reconstructed transcript and genes.bed is the reference. They look to be similar for the most part, but seems to have breaks between. What's the appropriate way of interpreting this result?

You've apparently colored your alignments by read strand. In this case, red indicates "+ (watson) strand" and blue indicates "- (strand) strand". This strand association is determined by the orientation of read #1 in a pair (or just "the read" if you have single-end data). This isn't actually documented in the user manual, but you can find one of the Broad institute software engineers replying here.

Your next question is likely why the color of the reads and the orientation of the genes are opposite of each other. This is because in modern directional protocols, read #2 in a pair matches the orientation of the gene from which it arose.


choosing a good representative genome subset I'm trying to build a genomic database for DNA alignments.
I started with NCBI accessions, but the data is very multiplicative, so I want to use subset of [max] N different strains for each specie.
my question is: I can i obtain the minimal N that the samples will be diverse enough and will cover the specie good?

Check out the progenomes database: http://progenomes.embl.de/