text:
Getting the licenses of all conda installed packages
I would like to get a list of the licenses of all packages I installed using Conda. pip has what I need in the form of pip-licenses, is there a similar thing for conda?

I proposed grep '"license":' ${CONDA_PREFIX}/conda-meta/*.json in a similar SO query.

To my knowledge there's no single listing of the licenses, but you can:
conda list | awk '{if(NR>3) printf("%s=%s", $1, $2)}' | xargs conda info | grep license

Allegedly conda info --license should work, but at least in conda 4.6.3 it doesn't show any license information. It unfortunate that there's no convenient way to get a tabular list of package:version:channel:license information from the command line.

It doesn't work for me, I get InvalidVersionSpec: Invalid version '0.1.0_license=1.1_r-mutex=1.0.0_tflow_180_select=3.0_tflow_select=2.1.0absl-py=0.2.2ad3=2.2.1alabaster=0.7.10anaconda-client=1.6.14anaconda-navigator=1.8.7anaconda-project=0.8.2anyqt=0.0.8appdirs=1.4.3argcomplete=1.9.4asn1crypto=0.24.0astor=0.6.2astroid=1.6.3astropy=3.0.2atomicwrites=1.1.5attrs=18.1.0awscli=1.15.57babel=2.5.3backcall=0.1.0backports=1.0backports.sh ... etc BTW, I can edit your response, isn't that strange? Conda version 4.6.3. Conda list does not seem to output any licenses.

Try adding a \n in the printf ("%s=%s\n"). Regarding the edit, it would become a suggestion and I'd get a notice for approval.

Thanx, this does produce a list of licences indeed! However, I don't see what they belong to. I also don't understand how your command works otherwise I would distill from it something that prints a packages name and a license only. What does the command do? There are no licences in the conda list output so where do they come from? How would I get the license for one package?

This one uses JQ for parsing and querying conda generated json output (https://stedolan.github.io/jq/)
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs -I{} conda search --info --json "{}" | jq --raw-output '.[][0] | "\(.name)\t\(.version)\t\(.license)"'

builds the list, asks conda search for all info on package as json, feeds json to JQ, which selects the first item and outputs its name, version and license as tsv.

I like this solution, although it does take a bit of time, since conda search queries aren't that fast. The only problem I ran into is if you have a later version of conda installed where the command line conda is actually a function, not an executable file, you'll need to changes ... xargs -I{} conda .. to ... xargs -I{} $CONDA_EXE ....

I can't comment, but here's an improvement for the top answer to list more information:
conda list | awk '{if(NR>3) printf("%s=%s\n", $1, $2)}' | xargs conda info | grep 'name\|version\|license'


text:
"Sequence Duplication Levels" module still fails after pre-processing Illumina data
I want to ask about why the sequence duplication levels are high after I trimmed by using Trimmomatic? I am using the following Trimmomatic operations: HEADCROP = 19 TRAILING = 20 MINLEN = 66.
How can i solve this problem? Thank You.

Why do you think this is a problem to begin with?

I though the cross sign (X) means some kind of error and should be eliminated with certain pre-processing technique (like trimming to solve adapter content error)? I am new to Illumina, thank you for advising :)

What type of sequencing libraries are these?  Whole genome, RNA Seq, whole exome, targeted seqeuencing?  What genome? also how complex is the library you sequencing.

To answer your direct question, there are a few reasons why there might be high levels of sequence duplication. From the FastQC help:

The underlying assumption of this module is of a diverse unenriched library. Any deviation from this assumption will naturally generate duplicates and can lead to warnings or errors from this module.

As @DevonRyan mentioned, with certain sequencing protocols such as RNA-Seq, two sequence reads at exactly the same location aren't that uncommon. This isn't a problem with RNA-Seq data, or with Trimmomatic, or with FastQC. It's just that this kind of data violates the assumption, and therefore should be ignored in those circumstances.
PCR duplicates are another possible cause. PCR duplicates can give the false impression of high coverage at a particular locus when in fact it's just a single observed read that has been duplicated many times (see here for more details). PCR duplicates can usually be detected and removed if your analysis involves mapping to a reference genome. But whether this is actually a problem you need to fix depends on what type of data you have and what types of analysis you want to do.
Large numbers of adapter dimers or rRNA may be present in your sample.

But I think it's also important to address how quality control (QC) is run. It can be tempting to run and re-run QC tools like Trimmomatic until all errors go away, but to be blunt these tools cannot think for you. For example, it's possible to get rid of most adapters by aggressively cropping/trimming both ends of each read, but you'll likely throw away a lot of good data that way. You may want to look into Trimmomatic's ILLUMINACLIP operation. It's also may be tempting to crop/trim reads aggressively if there are compositional biases near the beginning or end of the read. In fact, random hexamer priming can cause the Per-base Sequence Content module to fail on almost any RNA-Seq sample. Again, from the FastQC help (emphasis mine):

Libraries produced by priming using random hexamers (including nearly all RNA-Seq libraries) and those which were fragmented using transposases inherit an intrinsic bias in the positions at which reads start. This bias does not concern an absolute sequence, but instead provides enrichement of a number of different K-mers at the 5' end of the reads. Whilst this is a true technical bias, it isn't something which can be corrected by trimming and in most cases doesn't seem to adversely affect the downstream analysis. It will however produce a warning or error in this module.

So in other words, it's best to determine what can cause each FastQC module to fail, investigate whether this is actually a problem for your data set (referring to documentation as needed), and make a deliberate QC plan that addresses the issues that need attention.

FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether all of your samples are similar in their metrics.

source: Stack Exchange

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

sort: votes


text:
Why is Clostridium difficile toxin A + B outdated in snomed?
I'm looking at the new version of the Snomed dictionary and the concept 707993000 Clostridium difficile toxin A + B (substance) is deactivated (and outdated) since version 20170731. Can anybody explain me, why this happened?
We have used this concept for coding of results of laboratory examinations in our software system. Is there an alternative Snomed concept usable for the purpose of coding of evidence of the Clostridium difficile toxin A + B ?

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


text:
What is the purpose of the "fixedSearchFolder" and "fixedCombinedFolder" in MaxQuant?
What is the purpose of changing the fixedSearchFolder and fixedCombinedFolder in MaxQuant proteomics software? Is there any data that can be reused in other runs, so would it be good to use the same search folder for multiple runs? Or is it simply to direct temporary data and results data into different locations and each run should use a different folder?


text:
Meaning of the FORMAT fields of the VCF file coming from GIAB project
After reading the GIAB paper in https://www.biorxiv.org/content/early/2018/05/25/281006 and its Figure 1, I am still having trouble understanding the data inside the GIAB VCF file for HG001 (HG001_GRCh38_GIAB_highconf_CG-IllFB-IllGATKHC-Ion-10X-SOLID_CHROM1-X_v.3.3.2_all.vcf). In particular, I need to understand the information that led to a FILTER value of PASS for some variants.
The first case that I would like to discuss is the following variant:
chr1    110176467   .   A   G   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable    GT:PS:DP:ADALL:AD:GQ    0/1:.:512:114,106:145,148:760

Translating the information in the line to English, this is what I get:

We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
Across all the platforms, DP=512
ADALL (Values 114,106) is coming from “all" the datasets. Are those “all” the set (HiSeqPE300x + Cgnormal + 10XChromium), or the set (Illumina + CG + 10X + IonExome + SolidPE50x50bp + SolidSE75bp).
AD (Values 145,148) is coming from the same “all” datasets as ADALL, but they are unfiltered, so 145>114 and 148>106.

All this make sense, but I would like to know what case of Figure 1c was applied to arrive at the PASS value.
The second variant that I would like to mention is this one:
chr1    5705293 .   T   C   50  PASS    platforms=3;platformnames=Illumina,CG,10X;datasets=3;datasetnames=HiSeqPE300x,CGnormal,10XChromium;callsets=4;callsetnames=HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo;datasetsmissingcall=IonExome,SolidPE50x50bp,SolidSE75bp;callable=CS_HiSeqPE300xGATK_callable,CS_CGnormal_callable,CS_HiSeqPE300xfreebayes_callable;filt=CS_HiSeqPE300xGATK_filt   GT:PS:DP:ADALL:AD:GQ    0/1:.:599:121,139:49,28:627

Again, my translation would be:

We considered data from the following sequencing technologies: Illumina,CG,10X,IonExome,SolidPE50x50bp, and SolidSE75bp.
Of those, some of them (IonExome,SolidPE50x50bp,SolidSE75bp) were missing the call, so we discarded them.
Therefore, the data for the variant is going to come from 3 platforms (Illumina,CG,10X), 3 datasets (HiSeqPE300x,CGnormal,10XChromium), and 4 callsets (HiSeqPE300xGATK,CGnormal,HiSeqPE300xfreebayes,10XGATKhaplo), because we analyzed the data from HiSeqPE300x with GATK and also with Freebayes.
However, the callsets from CS_HiSeqPE300xGATK,CS_CGnormal, and CS_HiSeqPE300xfreebayes had this call in a region with low coverage of high MQ reads.
Also, the data from CS_HiSeqPE300xGATK_filt was filtered/discarded
Across all the platforms, DP=599
ADALL (Values 121,139) is coming from “all" the datasets. I assume (Cgnormal + 10XChromium), because CS_HiSeqPE300xGATK_filt was filtered.
AD (Values 49,28) are unfiltered.

But notice that 49<121 and 28<139. How come that there are less unfiltered reads (AD)? Should not be the case that AD > ADALL for all instances?
I don’t know how to reconcile the two variants. Aren’t they providing opposite information in the AD and ADALL fields? Which specific datasets and callsets are involved in getting AD and ADALL for each variant?
I would also appreciate answers to a couple of other questions:

Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?
Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

For reference, here are the descriptions in the header of the VCF file that I have been using:
##INFO=<ID=platforms,Number=1,Type=Integer,Description="Number of different platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=platformnames,Number=.,Type=String,Description="Names of platforms for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasets,Number=1,Type=Integer,Description="Number of different datasets for which at least one callset called this genotype, whether filtered or not">
##INFO=<ID=datasetnames,Number=.,Type=String,Description="Names of datasets for which at least one callset called this genotype, whether filtered or not">
##FORMAT=<ID=ADALL,Number=R,Type=Integer,Description="Net allele depths across all datasets">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Net allele depths across all unfiltered datasets with called genotype">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Total read depth summed across all datasets, excluding MQ0 reads">
##INFO=<ID=datasetsmissingcall,Number=.,Type=Integer,Description="Names of datasets that are missing a call or have an incorrect call at this location, and the high-confidence call is a variant">
##INFO=<ID=callable,Number=.,Type=String,Description="List of callsets that had this call in a region with low coverage of high MQ reads.">
##INFO=<ID=filt,Number=.,Type=String,Description="List of callsets that had this call filtered.">
##INFO=<ID=arbitrated,Number=1,Type=String,Description="TRUE if callsets had discordant calls so that arbitration was needed.">

The filtering criteria should be defined in the header of the VCF file. Can you include the header from your VCF file?
VCF Format and example header:
FILTERs that have been applied to the data should be described as follows:

##FILTER=<ID=ID,Description=”description”>

header:
##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
##INFO=<ID=AF,Number=.,Type=Float,Description="Allele Frequency">
##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
##INFO=<ID=DB,Number=0,Type=Flag,Description="dbSNP membership, build 129">
##INFO=<ID=H2,Number=0,Type=Flag,Description="HapMap2 membership">
##FILTER=<ID=q10,Description="Quality below 10">
##FILTER=<ID=s50,Description="Less than 50% of samples have data">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read Depth">
##FORMAT=<ID=HQ,Number=2,Type=Integer,Description="Haplotype Quality">

Thanks for the help. I posted the header fields that I am having trouble with.

For your second question, I don't think that variant matches any of the examples shown in Figure 1c. Note that that panel's caption says it shows four examples of arbitration. Per you last sub-question, I don't think this variant was subjected to arbitration.
For your second question, I believe the subtlety is in the field descriptions:
ADALL = Net allele depths across all datasets
AD = Net allele depths across all unfiltered datasets
with called genotype
My understanding is that ADALL is a count of allele depths that includes samples where the variant was filtered out or the genotype was not concordant; something to that effect. So generally I would expect ADALL > AD.
I'm honestly not sure exactly what they mean by "unfiltered datasets". That could mean "counting all reads without filtering on quality" or "callsets to which variant filtering was not yet applied". Maybe there's another subtlety there that explains why the two variants have opposite AD vs ADALL proportions.
For your additional questions:

Why are these two variants a PASS if there are 3 “callable” datasets in regions of low coverage?

According to the paper, having the variant in low-coverage regions isn't disqualifying if the genotype calls are concordant, which makes sense to me. Low coverage makes me more suspicious, but if multiple callsets are in agreement, I'm still willing to consider the call good (unless of course they all exhibit the same error mode, which is the main potential weakness of the consensus approach used by GiaB).

Did the variants follow an arbitration process? I guess not, if the INFO flag “arbitrated” is missing.

Right, if that flag is not present, it means arbitration was not needed/applied.

Thanks. I also think that there is a subtlety somewhere. Logically it does not make sense that ADALL >  AD, no matter if called or not.


text:
How to remove frame on RDKit figures?
RDKit can plot molecules, thought the structures are surrounded by an ugly frame with ticks. How can I turn that off and plot the molecules without that frame?
Note, the %pylab inline portion of this code is from the Jupyter environment.
%pylab inline
from rdkit import Chem
from rdkit.Chem import Draw

smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
title('Test')

It looks like rdkit and rdkit's Chem just use the matplotlib. So, if you import matplotlib you can edit your figure like any other matplotlib figure.
I produced this from the command line, so it may work differently in the Jupyter environment.
It looks like Chem doesn't properly scale the drawing based on how large the chemical is, so I added plt.xlim and plt.ylim to make the molecule larger in the frame.
from rdkit import Chem
from rdkit.Chem import Draw
import matplotlib.pyplot as plt
smiles = 'CCC'
m = Chem.MolFromSmiles(smiles)
fig = Draw.MolToMPL(m)
plt.title('Test')
plt.axis('off')
coord_min = 0.35
coord_max = 0.65
plt.xlim((coord_min, coord_max))
plt.ylim((coord_min, coord_max))
fig.savefig('mol.jpeg', bbox_inches='tight')


text:
What is the difference between a MeSH concept and a MeSH descriptor
Medical Subject Headings (MeSH) is a very comprehensive vocabulary to index the biomedical literature on Medline or PubMed (which is the major Biomedical literature database). MeSH is arranged into Descriptor > Concept > term hierarchy, but I still do not understand the difference between a MeSH descriptor and a MeSH concept.
For example, for Knee

MeSH descriptor is D007717
MeSH concept is M0012066 (Under the concepts section in the link)

Even after reading the description here and here, I am still not sure about the purpose of MeSH concept and MeSH descriptor.
Could anyone please shed some light on it?

MeSH is organized into Descriptor > Concept > Term hierarchy, where all the synonymous <term>'s are grouped together under a <concept>. For example, consider the MeSH terms,

Blood-Brain Barrier
Hemato-Encephalic Barrier
Brain-Blood Barrier

All these terms are group under the umbrella concept of "Blood-brain barrier".


text:
Machine learning using protein-sequences
I'm participating in a bioinformatics machine-learning seminar at my university. The main task is predicting binary classification of protein-protein interactions using sequence data as input.
One of the subtasks is familiarization with the dataset and presenting the dataset. Now I'm wondering which additional information I can get out of the sequence data.
I want to start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact. The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”. Species are completely mixed.
Is there something interesting that I can additionally include in my analysis? I'm used to exploratory data analysis though this is different due to the nature of the protein sequences.

Hi there, it would be perhaps good to describe the problem bit more. Do I get it right that you have a dataset of amino acid chains and now you want to predict witch do proteins interact? Do you have a gene position of proteins in the genome? Expression data? Do you have species? What is your training dataset?
These might be completely irrelevant questions since I am not a proteomic person. However, I don't really understand the question as it is right now.

Hi Kamil, thanks for the follow up question!
Yup, you get the problem right! I wanna start out with binary classification for the protein interaction with the labels “Interact” and “Non-Interact.
The dataset that I got provided consists of two fasta files. One containing the sequence header with the species and the corresponding amino-acid sequences ~300 entries. The other containing the exact same species header and the label “Interact” or "Non-interact”.
Species are completely mixed.

Great, we like juicy details. Could you perhaps [edit] your question and add it there?

Any specific kind of machine learning algorithm you intend or have to use?

Hi Florian, for our assignment we use ANN as well as SVM

There are two types of features which can be extracted from protein sequences.  It is not possible to known which, if any, of these features would be useful in classification tasks.  It may not be possible to build a classifier at all, or it may be very straight forward.  To know this, a feature selection technique must be used, such as e.g. forward feature selection.
Sequence based features
Sequence based features are derived from sequence-only information.
The most common used features are the PseAAC (pseudo amino acid composition) style features, including:

Average Position of each amino acid - 20 features
Average Distance between each amino acid - 20 features
Motifs - 20 features, 20*20 features, and/or 20*20*20 features (depending on motif length)
OAAC (Overall Amino Acid Composition) - 20 features
SAAC (Split Amino Acid Composition) - 60+ features
AAIndex Composition - 1 to 567+ features
PCC (Physicochemical Composition) - ~13 features
PSSM encoded as 1, 20, 200, 400 or 600 features (i.e. Psi-Blast output using NR or swissprot)
Average predicted residue burial/exposure level
Average predicted sequence entropy

Structure based features (predicted)
Using the sequence, you can either predict the 3d structure using homology, ab initio, or directly use a homologous structure.  You can then extract predicted 3d features, such as:

3d size (i.e. max of x,y,z)
Calculated ASA/SASA
3d contacts (type, frequency, etc)
any other local/global environmental information
cath/scop classification
dssp/stride composition


text:
structural variants: why doesn't AS + AP = SU?
I have run smoove then duphold on 46 individuals, then pasted the duphold output back into a squared vcf with all 46 samples.
I have been using SU/AC as one of my filters, looking for variants that have a minimum threshold of support per alternate allele called. Now I'm confused by some of the variants that are getting filtered out, and wonder why SU (number of pieces of evidence supporting the variant) is not calculated as I had expected.
For one variant, if I sum the AS values (split-read observation counts) and the AP values (paired-end observation counts) over all 46 individuals, I get AS=24 and AP=15, so 39 total reads supporting the variant in the population. However, SU is only 14 (PE=5 + SR=9).
Can anyone tell me why the values for SU, PE, SR are so much lower than the total AS & AP counts? Thank you!

We welcome this question here, but you may get a quicker and more thorough response posting to the smoove issue tracker.

Thanks - it's now posted here: https://github.com/brentp/smoove/issues/53


text:
Converting a BAM file into VCF
I have NGS illumina RNA-seq reads from M. musculus (mm10). I am trying to find variants along the strand portion of the reads in the refseq (mm10).
I mapped a pair of sequence files and generated a BAM file. Now I need to generate a VCF file from this BAM file. I’ve checked Galaxy for direct tools, but I couldn’t identify one. Any suggestions on how to turn a BAM file into VCF? Thank you in advance.

Possible duplicate of How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?

@gringer: The question asks about RNA-seq reads, where as your flagged potential duplicate is asking about DNA reads as best I can see. The process is not quite the same.

The approach works with RNASeq reads as well. All it needs is sufficient coverage across mapped reads. Where coverage is too high, reads can be sub-sampled via digital normalisation to reduce high-coverage regions without compromising coverage for low-coverage regions.

Freebayes giving me following error: terminate called after throwing an instance of 'std::out_of_range' what(): basic_string::substr: __pos (which is 166) > this->size() (which is 0)
Aborted (core dumped)

It's not really possible to convert bam to vcf. bam is a mapping file, it does not contain the information about variants, this information needs to be inferred in process called variant calling. I find important to mention that it's not just a different format of the same thing.
How to call variants (a vcf file) from mapped reads (a bam file) is very broad question, it depends on the sequencing technology (Illumina, PacBio, ...), the type of sequencing (whole genome sequencing, RAD, RNA-seq, ...), the sampling (calling of individuals vs population variability), the quality of a reference (human vs denovo assembly made by a postdoc in your lab)... For a lot of cases the answer is simply GATK, partially because it has well written manual, but mainly just because many people use it.
Some more detailed cases are well answered in more specific questions :

How do I generate a variant list (i.e. VCF file) using Illumina reads from a human genome?
Somatic tumor only variant calling?
Why are there missing calls in a VCF file from exome sequencing?

Indeed, relevant to this question, GATK have a Variants from RNA-seq best practices. Or they used to, and now it seems to have disappeared.

Honestly, I think good practices do not involve GATK. But it's true that GATK is what makes a review process much smoother (reviewers want GATK).

You can use Freebayes, provided you have your BAM file and the reference genome.
Example:
freebayes -f genome.fa aln.bam > var.vcf

I'd suggest you go through the whole github documentation, and choose the option(s) more suitable for your needs (ploidy, read support, etc.)

Thanks. Unfortunately when I tried this approach I didn’t get any variant readings along the mapped portion of the refseq. Not sure why. I know there are a few variants.

This answer would be more useful if it included a cautionary comment about how what the OP is asking about is not a format conversion, but an analytical process, ie variant calling, like Kamil S Jaron provided. In addition, the recommendation to "go through the whole github documentation" is unlikely to be satisfying given the OP is probably new to the field. Kamil's pointers to specific resources are much more helpful.

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


text:
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

Thanks for the nice solution. Is there any similar way to combine two vcf file for one subject (one vcf file is having SNP+Indel obtained from GATK HaplotypeCaller, GenotypeGVCF and another vcf file is having SVs from GenomeStrip). I'll definitely try this way as many researchers have suggested in Biostar also. Thanks.

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


text:
Merge 2 VCFs from different variant callers
I am working with WES data for detection of somatic variants and I have used two variant callers because no variant caller is complete in itself. I have used GATK Haplotypecaller for small variants like snp,ins,del etc. And for large variants I have used pindel (specially for SVs detection). Now I have two vcf file for each patients, one from gatk and another from pindel. Although I know various vcf-merging tools but is it right to merge two different vcf file from different callers. Do they give right and proper merged vcf file. Are those vcf-merging tools SV-aware and consider only trivial call for the same variant represented differently in different callers. I know these are too many questions for a single thread but all these are related to each other. I am in dilemma whether I should merge them or use them separately.
Thanks.

The merge utility from SURVIVOR seems to be what you are looking for. I did not try it yet, but it seems to be especially designed to handle SV.
The simplest thing to do would probably to first merge your patients files together for each tool and have one multi-patient file per variant caller. Then I would merge the multi-sample files of the different tools.
Here is a detailed blog post describing the different challenges associated with merging structural variants from different callers and explaining the steps they took to address it. All the code to reproduce the pipeline is available on a github repo.

Thanks for your reply and sorry for late response. Survivor is for merging SVs from different variant caller but can we merge small and large variants vcf file using this?

I am not sure, but It would seem it handles both small and large variants by looking at the source code. They also have a wiki detailing the methods. I think you might want to try on a test file and see if it does what you expect.

Thanks. I'll try and let you know if it works well with both variants or not.

Hey!! I know it is very late to respond but I was quite busy in other stuffs. I tried this tool and found that this tool extract SVs from both VCF file and merge them in a single file. But I need merging tool which combine all variants from both VCF file including SVs and merge than into single VCF file not just SVs. Unfortunately, that does not work. Any other idea... Thanks


text:
Protein sequence from patient data
Currently, I am working on NGS data and my aim is to get significance prediction of variants present in the vcf file. As we know about SIFT Score for significance score prediction, I am trying to understand how this score works. When I read its reference paper, in the method section, the procedure to determine the sift score starts with patient protein query sequence.
My question is how to get patient protein sequence from patient data (either vcf file or bam file) (which is not explained in the paper). Rest I can follow by other algorithms like PSI-BLAST & MOTIF.
Any suggestions....
Thanks

Ensembl now have a tool called Haplosaurus (still in beta) which will convert your phased VCF into actual protein sequences.

Haplosaurus
haplo is a local tool implementation of the same functionality that
  powers the Ensembl transcript haplotypes view. It takes phased
  genotypes from a VCF and constructs a pair of haplotype sequences for
  each overlapped transcript; these sequences are also translated into
  predicted protein haplotype sequences. Each variant haplotype sequence
  is aligned and compared to the reference, and an HGVS-like name is
  constructed representing its differences to the reference.
This approach offers an advantage over VEP's analysis, which treats
  each input variant independently. By considering the combined change
  contributed by all the variant alleles across a transcript, the
  compound effects the variants may have are correctly accounted for.
haplo shares much of the same command line functionality with vep, and
  can use VEP caches, Ensembl databases, GFF and GTF files as sources of
  transcript data; all vep command line flags relating to this
  functionality work the same with haplo.
Usage
Input data must be a VCF containing phased genotype data for at
  least one individual; no other formats are currently supported.
When using a VEP cache as the source of transcript annotation, the
  first time you run haplo with a particular cache it will spend some
  time scanning transcript locations in the cache.
./haplo -i input.vcf -o out.txt -cache

,@terdon : Thanks to both. I tried this. But I got some errors. Error is Attempt to reload Bio/EnsEMBL/VEP/TranscriptTree.pm aborted.. Although this file is present. I have checked  and confirmed this. I have also checked the similar issue at github (https://github.com/Ensembl/ensembl-vep/issues/169) which is still unsolved. Can you help me with this. Thanks.

You're best sticking with Irina on GitHub with this. She's working directly on developing the tool so will have the best insight and, if it's down to a bug, she can fix it.

I see that Irina has now solved this on GitHub. If anyone else has a similar problem, click through to the linked GitHub issue to see the resolution.

Thanks for software support. One last problem is stopping me to mark this as correct answer. I am not able to interpret output file. From your answer, I was expecting real protein sequence but I am getting many lines like this ENST00000342066 ENST00000342066:1027T>C ENSP00000342313:343W>R rs6672356 SM_5:2 in output file. Where is the protein sequence. How to read or interpret the output of ./haplo. Thanks.

The json option will give you sequence

Thanks for reply. I tried this too. The output file contains lots of sequences but it is very heavily packed with data or lots of data which is difficult to read. BTW, as you said, protein sequence can be obtained from phased vcf file which I got now. I have raised another issue related to performance of this tool please take a look at https://github.com/Ensembl/ensembl-vep/issues/214 . Thanks. Now I can mark it as right answer. Thanks again.

You can use SIFT in two different ways. One way is manually via their website. This is done by downloading the protein sequence (fasta) from a reliable repository (e.g., uniprot or NCBI), and change the amino acid of interest. Then calculate the SIFT score.
The other way, which is probably what you are interested in is with NGS data. Here normally SIFT is already known (precalculated) for all coding SNPs, and you can use that info directly during annotation of your variants. I don't know which variant annotation pipeline you are using, but for example annovar has SIFT data available at their databases repository.

You mean to say that for each variant at specific position (SNP especially), SIFT database has already determined which protein it will affect or what will be its corresponding protein sequence and how much similar it is with reference protein sequence and what will be its similarity score.  Is SIFT score valid for INDELs also.

Yes SIFT scores for the human genome are already precalculated, and you can use these scores for annotation of your variants (that is more efficient than to calculate them over and over again with each analysis, doesn't it?). SIFT works only for an amino acid substitution, that means nonsynonomous SNV. So no indels.


text:
Infer missing UTR features in GFF3 file
I am working on a GFF file that is missing the 5'UTR and 3'UTR information. For example:
  ctg123 . gene      1050  9000  .  +  .  ID=gene00001;Name=EDEN
  ctg123 . mRNA      1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
  ctg123 . exon      1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
  ctg123 . exon      3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
  ctg123 . exon      5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
  ctg123 . exon      7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
  ctg123 . CDS       1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1

Is there a way to add rows in the GFF with the 5'UTR and 3'UTR ranges?

Would we also need to deal with the negative strand? I mean, can you have gff lines referring to the negative strand and, if so, can we assume they will be correct? That the positions will be given with respect to the + strand and the end position will be "before" the start one?

yes, we do have to deal with the negative strand, too. Here is just a simplified example.

Then please [edit] your question and include an example of the negative strand as well. We would also need to know if each gff will have a single entry or if you can have multiple genes in the same gff. The more specific information you give us, the likelier it is that we will be able to give you a useful answer. However, bear in mind that UTRs can be spliced, so this is not a good way of getting the information you need. Getting a full annotation from wherever you downloaded this would be better.

Annoyingly, you can load such a GFF3 file into R (using the GenomicFeatures package) and access the UTRs, but they don't then get saved if you use the export function from rtracklayer (this would have been the simplest solution).
Since I don't know of anything to do this off-hand, I spent a couple minutes to write one. The usage would be addGFF3UTRs.py input.gff output.gff. You can modify a number of the labels to appear as you'd like. This hasn't been extensively tested, but it seems to produce the correct results in a small test I just ran.

The code is hosted on GitHub at the link I gave above and I also include it here:
#!/usr/bin/env python
import argparse

def parseAttributes(kvps):
    d = dict()
    for kvp in kvps:
        k, v = kvp.split("=")
        d[k] = v
    return d

def parseGFF3(fname, topLevel="gene", secondLevel="mRNA", CDS="CDS", exon="exon"):
    f = open(fname)
    genes = dict()
    transcripts = dict()
    g2t = dict()  # A dictionary with gene IDs associated to their transcript IDs

    for line in f:
        if line.startswith("#"):
            continue
        cols = line.strip().split("\t")
        attribs = parseAttributes(cols[8].split(";"))
        ID = attribs["ID"]

        if cols[2] == topLevel:
            genes[ID] = line
            g2t[ID] = []
        elif cols[2] == secondLevel:
            parent = attribs["Parent"]
            g2t[parent].append(ID)

            # Generate a list of [line, [exon entries], [CDS entries]]
            if ID not in transcripts:
                transcripts[ID] = [line, [], []]
            else:
                transcripts[ID][0] = line
        else:
            parent = attribs["Parent"]
            if parent not in transcripts:
                transcripts[parent] = [None, [], []]
            if cols[2] == exon:
                transcripts[parent][1].append((int(cols[3]), int(cols[4]), line))
            elif cols[2] == CDS:
                transcripts[parent][2].append((int(cols[3]), int(cols[4]), line))
    f.close()

    return genes, transcripts, g2t

def findUTRs(transcripts, fivePrime = True):
    for k, t in transcripts.items():
        # If there are no CDS entries then skip
        if len(t[2]) == 0:
            continue
        # Get the strand ("." will be treated as "+")
        strand = t[0].split("\t")[6]

        # For 3' UTR, just swap the strand
        if not fivePrime:
            strand = "+" if strand == "-" else "-"

        exons = [(s, e) for s, e, _ in t[1]]
        CDSs = [(s, e) for s, e, _ in t[2]]
        exons.sort()
        CDSs.sort()
        UTRs = []
        if strand != "-":
            final = CDSs[0][0]
            for s, e in exons:
                if e < final:
                    UTRs.append((s,e))
                elif s < final:
                    UTRs.append((s, final - 1))
                else:
                    break
        else:
            final = CDSs[-1][-1]
            for s, e in exons:
                if e < final:
                    continue
                elif s > final:
                    UTRs.append((s, e))
                else:
                    UTRs.append((final + 1, e))
        t.append(UTRs)

def saveGFF(oname, genes, transcripts, g2t, fivePrime, threePrime):
    o = open(oname, "w")
    o.write("##gff-version 3\n")
    for geneID, geneLine in genes.items():
        o.write(geneLine)
        # Go through each transcript
        for transID in g2t[geneID]:
            t = transcripts[transID]
            if t[0] is None:
                continue
            o.write(t[0])
            cols = t[0].strip().split("\t")

            # Write the exons, then CDS, then 5'UTR, then 3'UTR
            for exon in t[1]:
                o.write(exon[2])
            for CDS in t[2]:
                o.write(CDS[2])
            for idx, UTR in enumerate(t[3]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], fivePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.fivePrimeUTR{};Parent={}\n".format(transID, idx, transID))
            for idx, UTR in enumerate(t[4]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], threePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.threePrimeUTR{};Parent={}\n".format(transID, idx, transID))
    o.close()

parser = argparse.ArgumentParser(description="Parse a GFF3 file lacking UTR entries and add them in. Note that this program assumes a reasonably formatted GFF3 file")
parser.add_argument("--fiveUTRname", default="five_prime_UTR", help="The label for 5' UTR entries (default: %(default)s)")
parser.add_argument("--threeUTRname", default="three_prime_UTR", help="The label for 3' UTR entries (default: %(default)s)")
parser.add_argument("--topLevelID", default="gene", help="The 'type' designating the top-level entry (default: %(default)s)")
parser.add_argument("--secondLevelID", default="mRNA", help="The 'type' designating the second-level entry, typically something like 'mRNA' or 'transcript' (default: %(default)s)")
parser.add_argument("--exonID", default="exon", help="The 'type' designating exons (default: %(default)s)")
parser.add_argument("--CDSID", default="CDS", help="The 'type' designating CDS (default: %(default)s)")
parser.add_argument("input", metavar="input.gff3", help="Input file name")
parser.add_argument("output", metavar="output.gff3", help="Output file name")
args = parser.parse_args()

genes, transcripts, g2t = parseGFF3(args.input, topLevel=args.topLevelID, secondLevel=args.secondLevelID, exon=args.exonID, CDS=args.CDSID)

# Add the UTRs
findUTRs(transcripts)
findUTRs(transcripts, fivePrime=False)

# Write the output
saveGFF(args.output, genes, transcripts, g2t, args.fiveUTRname, args.threeUTRname)

Among other things, the canon-gff3 program in the AEGeAn Toolkit infers and prints missing features like UTRs.
$ canon-gff3 < genes.gff3 > genes-with-utrs.gff3

You can use the add_utrs_to_gff.py script provided by NCBI RefSeq for this purpose. Running the script on your example produces the following:
./add_utrs_to_gff.py example.gff3
ctg123  .  gene             1050  9000  .  +  .  ID=gene00001;Name=EDEN
ctg123  .  mRNA             1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
ctg123  .  exon             1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
ctg123  .  exon             3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
ctg123  .  exon             5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
ctg123  .  exon             7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
ctg123  .  CDS              1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  five_prime_UTR   1050  1200  .  +  .  ID=utr00002;Parent=mRNA00001
ctg123  .  three_prime_UTR  5301  5500  .  +  .  ID=utr00004;Parent=mRNA00001
ctg123  .  three_prime_UTR  7000  9000  .  +  .  ID=utr00005;Parent=mRNA00001


text:
In GFF3, annotating more than one protein-coding gene (i.e. polycistronic) contained in a eukaryotic mRNA
Below is an example of a simple GFF3 file:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS 3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS 3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

My question is: if we found another coding sequence (encoding a different protein) range from 3752 to 3904, how should the GFF3 be formatted? It seems to me that the GFF3 file can only allow one protein-coding gene per mRNA. If not, could anyone show me one example? Thank you!
[Edit]:
Sorry, using "CDS" is misleading. I changed the title and contents to "contain more than one protein-coding genes" instead of "more than one CDS". I also change the starting position of the second protein to 3752 (was 3750) to reflect the fact. My sincere apologies.

The GFF3 file format specification doesn't care where annotations are described in the file, although it can make some things easier if the annotations are ordered by the start position:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3750    3904    .   -   0   ID=ATRV01011:CDS:1;Parent=ATRV01011.0
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1

There are many examples in well-studied genomes of overlapping genes, so the GFF3 format needs to be flexible enough to handle these. Different genes can be encoded across the same region, but on different strands. But different genes can also be encoded on the same strand, sometimes sharing an exon.

thanks for your help. So essentially we can just use a new name (e.g here from ATNG to ATRV) to describe the new protein, correct? Also please see my edit. The second protein should actually be in a different open reading frame. Sorry I calculate it wrong the first time. However, I assume your answer still apply.

It is certainly possible to encode polycistrons in GFF3. In fact, this is discussed explicitly in the GFF3 specification.

Polycistronic transcripts
This is the case in which a single (possibly spliced) transcript encodes multiple open reading frames that generate independent protein products.
----->XXXXXXX*-->BBBBBB*--->ZZZZ*-->AAAAAA*-----
Since the single transcript corresponds to multiple genes that can be identified by genetic analysis, the recommended solution here is to create four "gene" objects and make them the parent for a single transcript. The transcript will contain a single exon (in the unspliced case) and four separate CDSs:
chrX  . gene XXXX YYYY  .  +  . ID=gene01;name=resA
chrX  . gene XXXX YYYY  .  +  . ID=gene02;name=resB
chrX  . gene XXXX YYYY  .  +  . ID=gene03;name=resX
chrX  . gene XXXX YYYY  .  +  . ID=gene04;name=resZ
chrX  . mRNA XXXX YYYY  .  +  . ID=tran01;Parent=gene01,gene02,gene03,gene04
chrX  . exon XXXX YYYY  .  +  . ID=exon00001;Parent=tran01
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene01
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene02
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene03
chrX  . CDS  XXXX YYYY  .  +  . Parent=tran01;Derives_from=gene04

To disambiguate the relationship between which genes encode which CDSs, you may use the Derives_from relationship.

This example doesn't indicate what the coordinates of each feature should be, but I think it makes most sense to have the gene and CDS feature coordinates correspond to the region coding each protein product, and the mRNA and exon feature coordinates correspond to the entire span.
BIG CAVEAT: Polycistrons are not typically considered a normal feature of eukaryotic genomes—they're much more common in bacteria and archaea. However, there are some confirmed cases of eukaryotic polycistronic mRNAs. That said, if overlapping gene predictions are observed in a GFF3 file, one cannot immediately assume that these are polycistronic. One or both of the gene predictions might be technical artifacts, or real genes that are incorrectly annotated, or (more rarely) actual overlapping genes. A claim that multiple protein products are expressed from a single transcript would require some convincing experimental validation.

If you found another CDS from 3750 to 3904 then it would necessarily be on a different transcript, given that you already have a CDS from 3860 to 3913 and having overlapping CDSs on the same transcript doesn't make much biological sense (though I'm sure there's some weird organism out there that has a gene like this). So the likely GFF3 file would look like:
1   T1  gene    3631    4605    .   +   .   ID=ATNG01010
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.1;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010:exon:1;Parent=ATNG01010.1
1   T1  CDS     3860    3913    .   +   0   ID=ATNG01010:CDS:1;Parent=ATNG01010.1
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010:exon:2;Parent=ATNG01010.1
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010:CDS:2;Parent=ATNG01010.1
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010:exon:3;Parent=ATNG01010.1
1   T1  mRNA    3631    4605    .   +   .   ID=ATNG01010.2;Parent=ATNG01010
1   T1  exon    3631    3913    .   +   .   ID=ATNG01010.2:exon:1;Parent=ATNG01010.2
1   T1  CDS     3759    3904    .   +   0   ID=ATNG01010.2:CDS:1;Parent=ATNG01010.2
1   T1  exon    3996    4276    .   +   .   ID=ATNG01010.2:exon:2;Parent=ATNG01010.2
1   T1  CDS     3996    4260    .   +   2   ID=ATNG01010.2:CDS:2;Parent=ATNG01010.2
1   T1  exon    4486    4605    .   +   .   ID=ATNG01010.2:exon:3;Parent=ATNG01010.2

Note that I've just added another transcript, since that's likely what would be biologically going on.

I calculate the ORF wrong (should be from 3752 to 3904, not 3750 to 3904). Please see the updated question. It does overlap with the original CDS but not in the same frame. This is a type of uORF. Any thought? Again, thanks for your help and sorry for the mistake.

Then the answer from @gringer is the more relevant one.

I found one similar example in an Arabidopsis annotation file. It is case of the cpuORF (CONSERVED PEPTIDE UPSTREAM OPEN READING FRAME ). Although the first and the third mRNA isoforms are identical, they repeated the mRNA isoform and name them differently to accommodate the additional ORF/CDS. I guess this is another way to do it.

See the ranges for AT5G03190.1 isoform and AT5G03190.3 isoform.
Chr5    TAIR10  gene    758374  760382  .   +   .   ID=AT5G03190;Note=protein_coding_gene;Name=AT5G03190
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.1;Parent=AT5G03190;Name=AT5G03190.1;Index=1
Chr5    TAIR10  protein 758793  760148  .   +   .   ID=AT5G03190.1-Protein;Name=AT5G03190.1;Derives_from=AT5G03190.1
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  five_prime_UTR  758374  758792  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  CDS 758793  760148  .   +   0   Parent=AT5G03190.1,AT5G03190.1-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.1
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.2;Parent=AT5G03190;Name=AT5G03190.2;Index=1
Chr5    TAIR10  protein 758539  760148  .   +   .   ID=AT5G03190.2-Protein;Name=AT5G03190.2;Derives_from=AT5G03190.2
Chr5    TAIR10  exon    758374  758660  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758539  758660  .   +   0   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  exon    758843  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  CDS 758843  760148  .   +   1   Parent=AT5G03190.2,AT5G03190.2-Protein;
Chr5    TAIR10  three_prime_UTR 760149  760382  .   +   .   Parent=AT5G03190.2
Chr5    TAIR10  mRNA    758374  760382  .   +   .   ID=AT5G03190.3;Parent=AT5G03190;Name=AT5G03190.3;Index=1
Chr5    TAIR10  protein 758539  758676  .   +   .   ID=AT5G03190.3-Protein;Name=AT5G03190.3;Derives_from=AT5G03190.3
Chr5    TAIR10  exon    758374  760382  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  five_prime_UTR  758374  758538  .   +   .   Parent=AT5G03190.3
Chr5    TAIR10  CDS 758539  758676  .   +   0   Parent=AT5G03190.3,AT5G03190.3-Protein;
Chr5    TAIR10  three_prime_UTR 758677  760382  .   +   .   Parent=AT5G03190.3


text:
Can a data file in VCF format be converted into FASTA?
I'm considering purchasing the 'MyGenome' product by Veritas Genetics to analyze my genome for a project. I'd like the data to be in FASTA format, but Veritas only provides VCF data. Is it possible to convert this VCF data into FASTA format?

What is "MyGenome" exactly? Will they sequence your entire genome? Only your exome? Specific regions of the exome?

@terdon Whole genome 30X on HiseqX. Link

Yeah, that doesn't sound like actual whole genome. Or, rather, it sounds like they might sequence the whole genome but then only analyze specific regions of it (which is a reasonable thing to do, I'm not saying they're swindling you). "myGenome is a whole genome screening test [. . . ]". If so, the method in b.nota's answer won't accurately reconstruct your genome, but only those regions of your genome that differ from the reference genome in the regions they happened t screen.

@terdon depends a bit on how they make the VCF file, don't you think? The interpretation is only based on specific regions, but they should have a file with all variants which OP should get for this VCF -> FASTA exercise.

@b.nota my guess (and it is only a guess) is that they will only do variant calling for targeted regions. That greatly speeds up the process and reduces the resources needed and also possibly protects them from litigation in some countries.

Good chance you guess right. OP might want to verify with MyGenome first, whether all variants are reported in VCF file.

Possible duplicate of How to manipulate a reference FASTA or bam to include variants from a VCF?

It's possible to create "a" fasta as discussed, but have my doubts whether that is a useful format for you to work with. But that's not the question, I guess.

You can try gatk function FastaAlternateReferenceMaker
java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \
   [--snpmask mask.vcf]

It might be worth pointing out that this will only reconstruct the OP's genome if the OP has WGS data. If not, I this tool will presumably use the reference genome for everything not explicitly mentioned in the file, so it won't be the OP's genome, per se.

@b.nota thank you for your answer. Just want to let you know that I've contacted veritas support for more information (still awaiting their response) and will accept your answer once I'm able to confirm that this works.


text:
CRISPR Sequence Finder and Database Download
I am searching for tools to pull CRISPR Spacers from Bacterial Genomes. I am aware of the CRISPRDB and the corresponding identification tool on the web server.

Are there other tools for finding CRISPR Spacers?
Is there a way to download the spacer sequences from CRISPRDB?
Is there another database that contains all validated / inferred bacterial CRISPR sequences?

CRT (CRISPR Recognition Tool) http://www.room220.com/crt/
Crisprs Finder Online Tool https://crispr.i2bc.paris-saclay.fr/Server/
Piler-CR http://www.drive5.com/pilercr/
(Omic Tools List) https://omictools.com/crispr-detection-category


text:
Where can I get the population allele frequency vcf file?
I want to use GATK to estimate cross-sample contamination for Whole Genome Sequencing data.
The specific tool is ContEst and it is run with:
java
 -jar GenomeAnalysisTK.jar \
 -T ContEst \
 -R reference.fasta \
 -I:eval tumor.bam \
 -I:genotype normal.bam \
 --popfile populationAlleleFrequencies.vcf \
 -L populationSites.interval_list
 [-L targets.interval_list] \
 -isr INTERSECTION \
 -o output.txt

The --popfile option requires a vcf file representing population allele frequencies. Where can I get it? Is it the one I can get from dbSNP?
I'm using the human_g1k_v37.fasta reference, a version of hg19, "for the phase1 analysis we mapped to GRCh37. Our fasta file which can be found here called human_g1k_v37.fasta.gz, it contains the autosomes, X, Y and MT but no haplotype sequence or EBV".

On the GATK forum they've recommended the population stratified VCF file for this purpose.


text:
What is the difference between a Bioinformatics pipeline and workflow?
I want to understand the difference between pipeline systems and workflow engines.
After reading A Review of Scalable Bioinformatics Pipelines I had a good overview of current bioinformatics pipelines. After some further research I found that there is collection of highly capable workflow engines. My question is then based on what I saw for argo. I would say I can be used as a bioinformatics pipeline as well.
So how do bioinformatics pipelines differ from workflow engines?

Great question! Note that from a prescriptive standpoint, the terms pipeline and workflow don't have any strict or precise definitions. But it's still useful to take a descriptive standpoint and discuss how the terms are commonly used in the bioinformatics community.
But before talking about pipelines and workflows, it's helpful to talk about programs and scripts. A program or script typically implements a single data analysis task (or set of related tasks). Some examples include the following.

FastQC, a program that checks NGS reads for common quality issues
Trimmomatic, a program for cleaning NGS reads
salmon, a program for estimating transcript abundance from NGS reads
a custom R script that uses DESeq2 to perform differential expression analysis

A pipeline or a workflow refers to a particular kind of program or script that is intended primarily to combine other independent programs or scripts. For example, I might want to write an RNA-seq workflow that executes Trimmomatic, FastQC, salmon, and the R script using a single command. This is particularly useful if I have to run the same command many times, or if the commands take a long time to run. It's very inconvenient when you have to babysit your computer and wait for step 3 to finish so that you can launch step 4!
So when does a program become a pipeline? Honestly, there are no strict rules. In some cases it's clear: the 10-line Python script I wrote to split Fasta files is definitely NOT a pipeline, but the 200-line Python script I wrote that does nothing but invoke 6 other bioinformatics programs definitely IS a pipeline. There are a lot of tools that fall in the middle: they may require running multiple steps in a certain order, or implement their own processing but also delegate processing to other tools. Usually nobody worries too much about whether it's "correct" to call a particular tool a pipeline.
Finally, a workflow engine is the software used to actually execute your pipeline/workflow. As mentioned above, general-purpose scripting languages like Bash, Python, or Perl can be used to implement workflows. But there are other languages that are designed specifically for managing workflows. Perhaps the earliest and most popular of these is GNU Make, which was originally intended to help engineers coordinate software compilation but can be used for just about any workflow. More recently there has been a proliferation of tools intended to replace GNU Make for numerous languages in a variety of contexts. The most popular in bioinformatics seems to be Snakemake, which provides a nice balance of simplicity (through shell commands), flexibility (through configuration), and power-user support (through Python scripting). Build scripts written for these tools (i.e., a Makefile or Snakefile) are often called pipelines or workflows, and the workflow engine is the software that executes the workflow.
The workflow engines you listed above (such as argo) can certainly be used to coordinate bioinformatics workflows. Honestly though, these are aimed more at the broader tech industry: they involve not just workflow execution but also hardware and infrastructure coordination, and would require a level of engineering expertise/support not commonly available in a bioinformatics setting. This could change, however, as bioinformatics becomes more of a "big data" endeavor.
As a final note, I'll mention a few more relevant technologies that I wasn't able to fit above.

Docker: managing a consistent software environment across multiple (potentially dozens or hundreds) of computers; Singularity is Docker's less popular step-sister
Common Workflow Language (CWL): a generic language for declaring how each step of a workflow is executed, what inputs it needs, what outputs it creates, and approximately what resources (RAM, storage, CPU threads, etc.) are required to run it; designed to write workflows that can be run on a variety of workflow engines
Dockstore: a registry of bioinformatics workflows (heavy emphasis on genomics) that includes a Docker container and a CWL specification for each workflow
toil: a production-grade workflow engine used primarily for bioinformatics workflows

Thanks, Daniel. Your answer outlines very well the details about the bioinformatics workflows. Your argument, that a very high level of software engineering is needed to use my listed workflow engines, is on point . If I understand correctly generally speaking bioinformatics pipelines can be pressed into workflow engines. So thanks again. That helped a lot.

Nice answer! I'm biased, but I feel like nextflow is also worth a mention :) And also WDL to complete the picture of common workflow tools. There's a nice "awesome-pipeline" repo list of workflow tools here too: https://github.com/pditommaso/awesome-pipeline - a little more comprehensive than the one in the original post, perhaps with more of a bioinformatics slant.

Thanks for the comment @tallphil. I too am biased by the tools with which I have experience. :-) Feel free to edit my answer to add a bit about nextflow and WDL.


text:
Gene function annotation - bacterial genome
I'm trying to annotate a genome to find all genes with a specific function. I have a FASTA and the read FASTQs - I'd like to assign the functional group of the identified proteins (e.g. Kegg orthology) automatically.
For more context, I have whole genome sequences of bacteria - I'm trying to find all genes involved in membrane transport. I'd either like a tool to retrieve gene ontologies from a list of gene names (input: genA output: Membrane Transport), or a tool that outputs the function as part of the annotation pipeline (e.g. genbank output with GO function defined).
Or posed in a different way, how would you identify all genes with a specific function from prokaryotic whole genome sequence data?
Are there any scripts that can do this?

Please [edit] your question and give us more details about what you have and what you need. Do you have gene models for the bacterium? I mean, do you first have to find the genes, then the gene functions, or do you already have the genes identified? You have FASTA and FASTQ of what? The entire genome? Raw reads? Are they assembled? Is there a reference sequence for the bacterium?

Thanks for your reply - edited in a few more details. I have the raw read data, and relative de novo assemblies of whole genome prokaryotes. I have a few species I'm looking at, a few E. coli, Staph, and Yersinia. Ref seqs available for all, but I'll be annotating draft clinical genomes.
I'm looking for a way to find all genes present in an assembly that match a specific function (e.g. membrane transport). In the genbank format for example I see "/function=" which is exactly what I want to assign - but I can't find a tool that writes that data. Prokka for e.g. does not write this by default

The PATRICBrc.org is a free resource you can annotate your genome and compare it with publicly available genomes.
Refer this tutorial for more detail.
Once you have annotated your genome in PATRIC, you will have PATRIC Local families (k-mer based protein family), which will help you identify all genes with a specific function.

Im not 100% on what you are looking for - you seem to jump from annotation to detection of genes in your question.
But if you have a database of genes you want to detect (presence/ absence) you could try srst2
which basically takes your raw reads (fastq) and maps them to the database.  it will give you a report on any hits, and any snps/ indels in the hit. The only problem is making you own database is difficult as you must adhere to the fasta header format.


text:
Calculating p-value for introns which are retained (or not) between 2 conditions
I am trying to calculate p-values and ultimately FDR values for RI introns for 2 conditions of experiment I have(control and knockdown)
For every condition I have 3 biological replicates each: 3 control replicates and 3 knockdown replicates.
This is the header line of my input file.
chrom   Istart  Iend    strand  avg-L_R_counts-1    Sk_counts-1 avg-L_R_counts-2    Sk_counts-2

-1 and -2 are 2 sample conditions control and knockout respectively
avg-L_R_counts = for a particular intron, this is the average of reads 5bp up stream and downstream of the 5' end and at the 3'end of the intron
Sk_counts = This is the number of reads spanning exon-exon junction for that intron.
I have approximately 350 introns(this is yeast). I can get the Percent Spliced In values but I am interested in getting a p-value and then I can calculate the FDR values.
Can anyone point me as to what approach I should take to calculate a p-value, or what test should I use?

What is your null hypothesis? What do you want to test? That one condition has more introns than the other?

My null hypothesis is that there is no intron retention in either conditions

Here is how we do this:
First break your gene model down into non-overlapping chunks, such that any possible gene model can be built from a combination of those chunks. E.g.
transcript 1 |>>>>>>>>>|----------|>>>>>>>>|---|>>>>|
transcript 2 |>>>>>>|-------------|>>>>>>>>|---|>>>>|
transcript 3 |>>>>>>|-------------|>>>>>>>>>>>>>>>>>|
chunks       |+++1++|+2|+++++3++++|+++4++++|+5+|+6++|

this can be achieved using cgat gtf2gtf --method=genes-to-unique-chunks from the cgat-apps toolkit if you don't want to code it yourself.
I call
chunks 1, 4 and 6 constitutive exon
chunk  3 constitutive intron
chunk  2 alternate intron
chunk  5 annotated retained intron.
Now use featureCounts to count the number of reads in each chunk.
Test the chunk counts for usage using DEXSeq.
Filter the results to just include chunks of type 2, 3 and 5 and recompute the FDRs.

If you have your transcripts in BED format, you can use bedops --partition to make chunks.

Did you look into literature? I am not an expert in this field, but a small search on google leads me to IRFinder (ref.), a tool that can assess intron retention. My advice is to see if this tool can help you, before reinventing the wheel...

You might take a look at the MISO algorithm for inspiration, or as a method for calculating confidence intervals of alternatively spliced isoforms: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3037023/


text:
Definition of "seed" in sequence alignment
I would like to know what is meant by "seed" for various sequence aligners. How is it important?

The seed is the subset of a read used in the first step of an alignment. Many aligners work by a seed-and-extend model, wherein they first find all regions matching the "seed" and then extend the alignment around that allowing mistmatches and indels until it either gives up (and therefore uses a different seed) or finds a sufficiently good alignment.

Thanks devon for the explanation. I have a slightly similar question as to why finding a seed for a very short exons is difficult? Considering a 3 nt exon, is it because seed  such small has a probability of mapping every 4^3=64 bases and virtually impossible to look for in the genome? Your views on this

Usually seeds are a few times that length, since you want few enough hits to limit your search space quickly.

Devon's answer gives a good, concise definition. But it's also helpful to consider why seed-and-extend is used and what benefits it provides.
Finding approximate string matches requires operations that are expensive computationally. On the other hand, finding exact string matches can be performed using much less expensive operations—that is, it can be done more quickly with less memory. So quickly finding the location of all seed matches and then applying the approximate matching ONLY AT THOSE LOCATIONS (the seed-and-extend strategy) makes sure that less time is wasted on chasing spurious similarities.
So how long should a seed be? Riffing off your comment, a seed 3bp in length would be entirely ineffective since that same 3bp sequence is going to occur ALL OVER THE GENOME IN RANDOM PLACES: genes, pseudogenes, transposable elements, regulatory sequences, you name it. On the other hand, a seed 51bp in length would also be ineffective since almost every seed would contain at least one sequencing error or variant/polymorphism, thus preventing a perfect match.
In practice, mappers designed for aligning 100-300bp Illumina reads typically use seeds in the high teens or low 20s.


text:
Find overlap between VCF files
I have two VCF files and I want to compare the missing rate in each of them, but I want to only look at sites that are present in both files.  How could I go about getting a list of positions that are included in both VCF files and then filtering those VCF files to only include those positions?

Something like vcf-isec? Quoting:

Creates intersections and complements of two or more VCF files. Given multiple VCF files, it can output the list of positions which are shared by at least N files, at most N files, exactly N files, etc. The first example below outputs positions shared by at least two files and the second outputs positions present in the files A but absent from files B and C.

Thank you!!!  I had suspected that there would be something like that but I was having trouble finding anything. I really appreciate it

The bedtools intersect command is my go-to for operations like this. From the manual:

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

You can try something like this in bcftools, which is a bit quicker than VCFtools, as it uses block gzipping and indexing.
# make sure both files are indexed and in bgzip format
bcftools index vcfA.vcf.gz
bcftools index vcfB.vcf.gz

# find overlapping positions between the two files
bcftools isec -n~11 -c all vcfA.vcf.gz vcfB.vcf.gz > sharedPositions.txt

# filter vcf files for the overlapping positions
bcftools view -T sharedPositions.txt vcfA.vcf.gz -Oz > vcfA.sharedPositions.vcf.gz

bcftools view -T sharedPositions.txt vcfB.vcf.gz -Oz > vcfB.sharedPositions.vcf.gz

The perhaps use something like plink2 to look at the missing rate per variant and individual for each file:
plink2 --vcf vcfA.vcf.gz --missing --out vcfA.missingStats


text:
Any fast options to query large VCF bed intervals?
I'm doing some analysis and I need to subset a large VCF file (~8GB gziped) given a bed interval and identify within a subset of rsid.
Unfortunately, both my normal choices to do this analysis (snpSift and bedtools) are taking way to long or failing due to memory issues in my local computer and a remote server.
Do you guys know any other options or suggestions to speed up this process?
Follow the commands I use:
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

or
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

The bedtools command usually fail due to unknown reason and SnpSift runs over 6 hours even given 10GB of ram. My local machine have 8GB of RAM, but the server have 32GB.

I forgot to ask about versions. Maybe a newer version is faster ? What are those "unknown reasons" is there any error or message or just a freezing computer? (Maybe there's a bug). In the second command can't you combine the two snpSift commands in one? How big is the region you want?

@Llopis I'm using the latest versions of both programs. The bedtools error is usually unrecognized format and the command fails.
As far as I know, I can't combine the commands but I never had problem with this scheme.
And the region was approximately 32-40Mb

this:
bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]"

can be replaced with one GATK SelectVariants https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php
java -jar gatk.jar -T SelectVariants -R ref.fa \
  -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf

and furthermore, you can also work in parallel by splitting your bed file (e.g. by chromosomes), and merge the results.

thanks for the suggestion, but on my experience GATK usually takes longer to run this kind of process.
I've also been able to complete the filtering by inverting the snpSift selection order.

@andre if you could post the solution it might help other people

@Llopis sure, please find it bellow

Generically, with BEDOPS:
$ vcf2bed < <(gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed

Or:
$ vcf2bed < <(gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\t' myRegions.bed - > answer.bed

Etc.
BEDOPS supports input streams; use them where you can.
That said, if you're going to query the VCF data frequently, convert once and use it in downstream operations.
$ vcf2bed < <(gunzip -c snps.vcf) > snps.bed

Then:
$ bedops -e 1 snps.bed myRegions.bed > answer.bed

Etc.
If you have a computational cluster, you can further parallelize work by splitting SNPs on chromosomes with bedextract:
$ for chr in `bedextract --list-chr snps.bed`; do echo ${chr}; bedextract ${chr} snps.bed > snps.${chr}.bed; done

Then submit a batch job for each chromosome:
$ bedmap --chrom ${chr} --echo --echo-map-id --delim '\t' myRegions.bed snps.${chr}.bed > answer.${chr}.bed

Using bedextract and adding --chrom ${chr} to your bedmap call will yield large speed improvements over other toolkits, as these options focus work on the chromosome of interest and nothing else.

I was tinkering with the command and was able to complete the execution inverting the order of rsid and bed intervals filtering. The command is as follows:
gzcat <myvcf>.vcf.gz | \
    java -Xmx10g -jar snpSift.jar filter --set <myrsid>.txt "ID in SET[0]" | \
    java -Xmx10g -jar snpSift.jar intervals <mybed>.bed

Maybe my rsid subset is smaller than my BED intervals and thus the process is quicker and more efficient.
Anyway thanks everyone

Here is another way to do it in 2019:
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf

Things can be sped up if we compress and index the input.vcf. Then bcftools will have near-immediate random access to the position.
$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz
$ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf.gz > output.vcf

Internally bcftools works with bcf. Whenever the input and output is vcf it converts vcf>bcf>vcf. When dealing with large files, this produces a large overhead. Converting vcf to bcf first can reduce runtime for further downstream analyses dramatically. In my experience 2-3x faster.
$ bcftools view -Ob input.vcf > input.bcf
$ bcftools index input.bcf
$ bcftools view -Ob -T regions.bed -i 'ID=@<myrsid>.txt' input.bcf > output.bcf

Helpful comment about VCF > BCF > VCF conversion!


text:
Retrieving a list of human genes having GO associations
I need a command in R to retrieve all human genes associated with a Gene Ontology entry. I tried to look for it online but did not find it.

@DanielStandage I am not sure this question is about the human genome.

@Llopis When he says "retrieve all human genes" it definitely makes me think of the human genome. :-)

@Daniel yes, but the question is not about the human genome. Or should we also use gene as tag here ? Maybe we use different definitions of the tag

Here's an example for the mouse genome:
library(org.Mm.eg.db)
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))

You get output like:
                                                    GENENAME SYMBOL
1                                     pregnancy zone protein    Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
3 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
4                                             neurotrophin 3   Ntf3
5             neurotrophic tyrosine kinase, receptor, type 1  Ntrk1
6            furin (paired basic amino acid cleaving enzyme)  Furin
7              proprotein convertase subtilisin/kexin type 6  Pcsk6
8                                                 sortilin 1  Sort1
9                                      alpha-2-macroglobulin    A2m

Edit: Given your comment that you want instead a list of all genes with a GO term associated with them (for human):
library(org.Hs.eg.db) # Install it from bioconductor yourself
unlist(as.list(org.Hs.egSYMBOL)[mappedkeys(org.Hs.egGO)])

The mappedkey() part gets the index into the bimap of all entries with a valid GO mapping. That's then used to subset the symbol bimap (after converting to a list).

I guess the following code will help,
source("https://bioconductor.org/biocLite.R")
biocLite("biomaRt")
library("biomaRt")
ensembl = useMart("ensembl", dataset = "hsapiens_gene_ensembl")
#listAttributes(ensembl)
mapping <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol", "go_id"), mart = ensembl)
head(mapping)


text:
Error "Can't locate List/MoreUtils.pm in @INC" when I run update_blastdb.pl
I encoutered some technical problem when learning BLAST as below:
After I use conda install blast to install ncbi-blast, and run the command:
update_blastdb.pl --decompress 16SMicrobial
the shell complains:

Can't locate List/MoreUtils.pm in @INC (you may need to install the
  List::MoreUtils module) (@INC contains: /etc/perl /usr/loc
  al/lib/x86_64-linux-gnu/perl/5.22.1 /usr/local/share/perl/5.22.1
  /usr/lib/x86_64-linux-gnu/perl5/5.22 /usr/share/perl5 /usr/li
  b/x86_64-linux-gnu/perl/5.22 /usr/share/perl/5.22
  /usr/local/lib/site_perl /usr/lib/x86_64-linux-gnu/perl-base .) at
  /home/mck f111/miniconda3/envs/bioinfo/bin/update_blastdb.pl line 43.

After I run this command:
cpanm List::MoreUtils

It appears some installation happened, however, the complains are still there.
Could anyone tell me how to fix this? Thanks a lot!

What operating system are you running this on? That looks like some sort of *nix, but is it Linux, macOS, UNIX, Ubuntu subsystem on Windows? Something else? Did you run the cpanm command as a regular user or as root? Did you see where it installed (it should be mentioned in the command's output)? Please [edit] your question and clarify.

Install List::MoreUtils using the same user as blast.


text:
Find intervals that genes fall within their range
I am trying to run some frequency based stats to identify selective sweep in my study system using Rpackage "PopGenome".
To proceed with it I have split my genome data (whole genome sequencing) into 5MB chunks
cregions <- character(31)
cregions[1] <- "1-5000000"
cregions[2] <- "5000001-10000000"
cregions[3] <- "10000001-15000000"
cregions[4] <- "15000001-20000000"
cregions[5] <- "20000001-25000000"
cregions[6] <- "25000001-30000000"
cregions[7] <- "30000001-35000000"
....

I have a bunch of genes and I want to know each of my genes falls within the range of which intervals in the above list.
these are my genes:
chrom start stop gene_id
Chr01 107310787 107323633 gene2748
Chr01 113789595 113794851 gene2894
Chr01 159041389 159068535 gene4806
Chr02 2742629 2758060 gene4846

I'd like to get an out put like this
gene2748  window_2
gene2894  window_3
...

Any idea how I can do it in R?

The GenomicRanges package is a great place to start, something like this:
library(GenomicRanges)

# make a GenomicRanges set corresponding to the genes
genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom")

# turn your windows into a data frame and make a GenomicRanges set from that
windows = ...
windows.gr = makeGRangesFromDataFrame(df = windows ...

# find overlaps between the two sets
olaps = findOverlaps(genes.gr, windows.gr)

# output relevant fields
cat(paste(genes[queryHits(olaps),"gene_id"], windows[subjectHits(olaps),"window_id"]))


text:
Bash script error at paste command
I wrote script for pasting rsids on CADD output. Here is script.
#!/bin/bash
cd tmp
cut -f 1,2 CADD.tsv > fileA
paste fileA <(cut -f 2,125  CADD.tsv) > myNewFile
bedtools intersect -a  myNewFile -b New.vcf -wb |cut -f 1-4,7 > CADD.rsids.tsv

I have tested commands one by one and they work fine. But when I run script of these commands, it gives me this error.
$ sh cadd.rsids.sh
4: cadd.rsids.sh: Syntax error: "(" unexpected

Can anyone please tell me how to resolve this error?
Here are few lines of CADD.tsv file.
"##CADD GRCh38-v1.4 (c) University of Washington, Hudson-Alpha Institute for Biotechnology and Berlin Institute of Health 2013-2018. All rights reserved."
#Chrom  Pos Ref Alt Type    Length  AnnoType    Consequence ConsScore   ConsDetail  GC  CpG motifECount motifEName  motifEHIPos motifEScoreChng oAA nAA GeneID  FeatureID   GeneName    CCDS    Intron  Exon    cDNApos relcDNApos  CDSpos  relCDSpos   protPos relProtPos  Domain  Dst2Splice  Dst2SplType minDistTSS  minDistTSE  SIFTcat SIFTval PolyPhenCat PolyPhenVal priPhCons   mamPhCons   verPhCons   priPhyloP   mamPhyloP   verPhyloP   bStatistic  targetScan  mirSVR-Score    mirSVR-E    mirSVR-Aln  cHmm_E1 cHmm_E2 cHmm_E3 cHmm_E4 cHmm_E5 cHmm_E6 cHmm_E7 cHmm_E8 cHmm_E9 cHmm_E10    cHmm_E11    cHmm_E12    cHmm_E13    cHmm_E14    cHmm_E15    cHmm_E16    cHmm_E17    cHmm_E18    cHmm_E19    cHmm_E20    cHmm_E21    cHmm_E22    cHmm_E23    cHmm_E24    cHmm_E25    GerpRS  GerpRSpval  GerpN   GerpS   tOverlapMotifs  motifDist   EncodeH3K4me1-sum   EncodeH3K4me1-max   EncodeH3K4me2-sum   EncodeH3K4me2-max   EncodeH3K4me3-sum   EncodeH3K4me3-max   EncodeH3K9ac-sum    EncodeH3K9ac-max    EncodeH3K9me3-sum   EncodeH3K9me3-max   EncodeH3K27ac-sum   EncodeH3K27ac-max   EncodeH3K27me3-sum  EncodeH3K27me3-max  EncodeH3K36me3-sum  EncodeH3K36me3-max  EncodeH3K79me2-sum  EncodeH3K79me2-max  EncodeH4K20me1-sum  EncodeH4K20me1-max  EncodeH2AFZ-sum EncodeH2AFZ-max EncodeDNase-sum EncodeDNase-max EncodetotalRNA-sum  EncodetotalRNA-max  Grantham    Dist2Mutation   Freq100bp   Rare100bp   Sngl100bp   Freq1000bp  Rare1000bp  Sngl1000bp  Freq10000bp Rare10000bp Sngl10000bp EnsembleRegulatoryFeature   dbscSNV-ada_score   dbscSNV-rf_score    RemapOverlapTF  RemapOverlapCL  RawScore    PHRED
1   3362704 T   A   SNV 0   Transcript  INTRONIC    2   intron  0.536423841 0.04    NA  NA  NA  NA  NA  NA  ENSG00000142611 ENST00000270722 PRDM16  CCDS41236.2 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  118572  24215   NA  NA  NA  NA  0   0   0   -0.405  -1.689  -1.593  964 NA  NA  NA  NA  0   1   1   5   0   0   1   0   0   0   0   1   0   1   12  0   0   0   0   0   8   11  5   2   0   NA  NA  2.46    -4.93   1   0.02    4.53362 1.5498  3.49082 1.00224 9.30658 3.84201 5.88121 2.95189 9.31556 5.27366 4.01294 1.08327 25.99   5.10373 4.47027 1.51567 1.55624 0.7046  4.95028 1.49632 8.57929 1.45764 1.19689 0.644106    0.01783 0.01476 NA  16  1   0   15  4   10  184 17  101 1741    NA  NA  NA  5   19  -0.173019   0.424
1   7785635 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049245 ENST00000054666 VAMP3   CCDS88.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7785635 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.483443709 0.053333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  112 2704    NA  NA  NA  NA  0.009   0   0   0.418   -0.242  -0.522  726 NA  NA  NA  NA  0   0   0   3   1   0   14  1   0   0   1   0   0   0   22  0   0   0   0   0   1   1   4   0   0   NA  NA  5.52    -11 NA  NA  20.4013 3.22673 13.621  2.38936 11.7271 2.25561 9.94895 2.15059 16.5529 2.79938 11.0724 2.29032 12.2962 2.07181 32.6588 8.63119 15.114  2.76883 7.39549 2.14311 14.1811 3.08082 0.603948    0.246907    0.55734 0.16905 NA  9   1   1   15  3   5   174 35  68  1667    NA  NA  NA  2   4   0.049024    1.923
1   7803233 C   T   SNV 0   Transcript  INTRONIC    2   intron  0.331125828 0.013333333 NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  182 726 NA  NA  NA  NA  0.014   0   0   -0.553  -1.498  -0.938  673 NA  NA  NA  NA  0   0   0   6   0   0   23  5   0   0   0   0   0   0   14  0   0   0   0   0   0   0   0   0   0   NA  NA  7.49    -15 NA  NA  4.36329 1.02663 4.55863 1.4118  4.72226 1.51608 6.80823 2.09997 7.44497 2.86816 5.71457 2.35876 3.06623 0.62198 18.3251 3.23263 2.1423  0.84413 8.63178 2.70111 6.29127 1.59229 0.198671    0.048037    0.0465  0.01705 NA  6   1   1   15  2   2   134 26  39  1314    NA  NA  NA  NA  NA  -0.395046   0.067
1   7808665 T   C   SNV 0   Intergenic  DOWNSTREAM  1   downstream  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7808665 T   C   SNV 0   Transcript  INTRONIC    2   intron  0.344370861 0   NA  NA  NA  NA  NA  NA  ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 12:00:00 AM NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  5614    1683    NA  NA  NA  NA  0.082   0   0   0.463   -0.365  -0.373  659 NA  NA  NA  NA  0   0   0   5   3   1   24  0   0   0   0   0   0   0   15  0   0   0   0   0   0   0   0   0   0   NA  NA  5.23    -10.5   NA  NA  7.16527 1.57568 5.63104 1.29085 5.96112 1.41329 8.31159 2.11434 4.41828 1.14843 5.03919 1.48679 3.91787 1.13788 4.97324 1.34708 5.22225 1.10759 7.38244 1.40721 11.4505 2.73652 0.977252    0.19269 0.09436 0.03409 NA  43  1   0   10  3   2   110 30  49  1288    NA  NA  NA  7   26  0.11875 2.796
1   7827519 C   G   SNV 0   CodingTranscript    NON_SYNONYMOUS  7   missense    0.529801325 0.066666667 NA  NA  NA  NA  P   A   ENSG00000049246 ENST00000613533 PER3    CCDS72695.1 NA  18/22   2854    0.45172523  2590    0.712909441 864 0.714049587 lcompl  NA  NA  177 17277   tolerated   0.05    benign  0.127   0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 27  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7827519 C   G   SNV 0   RegulatoryFeature   REGULATORY  4   regulatory  0.529801325 0.066666667 NA  NA  NA  NA  NA  NA  NA  ENSR00000000832 NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  177 17277   NA  NA  NA  NA  0.041   0   0   -0.553  -0.852  0.026   525 NA  NA  NA  NA  1   0   0   3   0   2   8   9   1   1   1   0   0   1   2   6   1   2   9   0   0   0   0   1   0   NA  NA  12.1    -24.2   NA  NA  8.89687 3.01407 58.9843 24.0334 58.3619 18.9077 55.4967 17.2149 10.5523 3.06178 87.6613 48.8498 8.66388 2.84705 31.4762 9.05454 11.8069 3.99886 8.7429  2.26128 72.4647 18.4371 1.08886 0.638979    2.10021 1.18017 NA  2   2   1   19  4   8   204 31  59  1420    Promoter    NA  NA  86  94  0.346155    6.007
1   7828155 A   G   SNV 0   Intergenic  UPSTREAM    1   upstream    0.364238411 0.066666667 NA  NA  NA  NA  NA  NA  ENSG00000236266 ENST00000451646 Z98884.1    NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  NA  813 17022   NA  NA  NA  NA  0.009   0   0   -0.389  -0.735  -0.731  522 NA  NA  NA  NA  0   0   0   4   0   0   20  4   1   0   0   0   0   0   14  1   0   0   3   0   0   0   0   1   0   NA  NA  4.56    4.56    1   0.12    10.0753 2.35506 38.5512 20.107  21.923  10.3608 21.217  11.8934 5.61317 1.38406 9.74608 6.16783 3.76905 1.10216 24.6436 5.18488 14.8421 7.93141 4.56698 1.13989 18.7058 4.31099 0.500167    0.211141    0.16854 0.06818 NA  20  2   0   14  4   4   159 34  57  1421    NA  NA  NA  11  14  0.225659    4.402

Here are lines of New.vcf
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO
1   3362704 rs11807862  T   A   923.01  PASS    "BaseCounts=31,0,0,30;BaseQRankSum=-0.108;DB;Dels=0;FS=3.561;GC=61.35;HaplotypeScore=1.7256;MQ=60;MQ0=0;MQRankSum=-0.368;QD=15.13;ReadPosRankSum=0.498;DP=282;AF=0.5;MLEAC=1;MLEAF=0.5;AN=14;AC=7"
1   7785635 rs228729    T   C   2294.01 PASS    "BaseCounts=0,44,0,0;DB;Dels=0;FS=0;GC=42.39;HaplotypeScore=0;MQ=60;MQ0=0;QD=28.13;BaseQRankSum=2.114;MQRankSum=-1.268;ReadPosRankSum=-0.548;DP=1430;AF=0.5;MLEAC=1;MLEAF=0.5;AN=70;AC=51"
1   7803233 rs228642    C   T   2082.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=-0.963;DB;Dels=0;FS=2.884;GC=44.89;HaplotypeScore=0.9999;MQ=60;MQ0=0;MQRankSum=0.058;QD=14.92;ReadPosRankSum=0.409;DP=1468;AF=0.5;MLEAC=1;MLEAF=0.5;AN=62;AC=44"
1   7808665 rs228666    T   C   1925.01 PASS    "BaseCounts=0,36,0,26;BaseQRankSum=0.064;DB;Dels=0;FS=0.979;GC=30.42;HaplotypeScore=1.7333;MQ=60;MQ0=0;MQRankSum=-0.82;QD=18.03;ReadPosRankSum=-1.32;DP=939;AF=0.5;MLEAC=1;MLEAF=0.5;AN=42;AC=25"
1   7827519 rs228697    C   G   786.01  PASS    "BaseCounts=0,14,24,0;BaseQRankSum=0.016;DB;Dels=0;FS=4.925;GC=58.35;HaplotypeScore=0;MQ=60;MQ0=0;MQRankSum=-0.698;QD=18.63;ReadPosRankSum=0.99;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=9"
1   7828155 rs2859388   A   G   2385.01 PASS    "BaseCounts=0,0,56,0;DB;Dels=0;FS=0;GC=35.41;HaplotypeScore=0;MQ=60;MQ0=0;QD=24.91;BaseQRankSum=-0.231;MQRankSum=0.948;ReadPosRankSum=-0.128;DP=1136;AF=0.5;MLEAC=1;MLEAF=0.5;AN=46;AC=31"
1   7830057 rs2640909   T   C   1421.01 PASS    "BaseCounts=0,19,0,19;BaseQRankSum=1.051;DB;Dels=0;FS=5.041;GC=56.11;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=0.496;QD=16.61;ReadPosRankSum=-0.526;DP=638;AF=0.5;MLEAC=1;MLEAF=0.5;AN=32;AC=19"
1   11107089    rs12139042  G   A   1668.01 PASS    "BaseCounts=43,0,0,0;DB;Dels=0;FS=0;GC=45.64;HaplotypeScore=0.734;MQ=60;MQ0=0;QD=29.7;BaseQRankSum=-1.652;MQRankSum=0.027;ReadPosRankSum=-1.192;DP=354;AF=0.5;MLEAC=1;MLEAF=0.5;AN=18;AC=10"
1   11787392    rs3737967   G   A   719.01  PASS    "BaseCounts=23,0,23,0;BaseQRankSum=-0.475;DB;Dels=0;FS=4.566;GC=61.85;HaplotypeScore=0.734;MQ=60;MQ0=0;MQRankSum=-0.84;QD=15.63;ReadPosRankSum=0.329;DP=378;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11790870    rs2274976   C   T   777.01  PASS    "BaseCounts=0,18,0,24;BaseQRankSum=-0.042;DB;Dels=0;FS=4.677;GC=56.11;HaplotypeScore=0.6651;MQ=60;MQ0=0;MQRankSum=1.204;QD=15.05;ReadPosRankSum=0.836;DP=355;AF=0.5;MLEAC=1;MLEAF=0.5;AN=20;AC=10"
1   11792243    rs1476413   C   T   1751.01 PASS    "BaseCounts=1,18,0,21;BaseQRankSum=1.028;DB;Dels=0;FS=0;GC=56.11;HaplotypeScore=1.7287;MQ=60;MQ0=0;MQRankSum=-0.408;QD=16.65;ReadPosRankSum=-0.099;DP=954;AF=0.5;MLEAC=1;MLEAF=0.5;AN=50;AC=29"
1   11794400    rs4846051   G   A   2616.01 PASS    "BaseCounts=46,0,0,0;DB;Dels=0;FS=0;GC=59.35;HaplotypeScore=6.2573;MQ=60;MQ0=0;QD=31.75;DP=1654;AF=1;MLEAC=2;MLEAF=1;AN=78;AC=78"
1   11794419    rs1801131   T   G   2028.01 PASS    "BaseCounts=0,1,21,21;BaseQRankSum=-0.117;DB;Dels=0;FS=9.995;GC=59.35;HaplotypeScore=2.4328;MQ=60;MQ0=0;MQRankSum=1.03;QD=14.84;ReadPosRankSum=0.665;DP=1176;AF=0.5;MLEAC=1;MLEAF=0.5;AN=54;AC=33"

The error complains about the parenthesis "("

If I remove "(" parathesis, this gives me error of 4: c.rsids.sh: cannot open cut: No such file

Try to write a pipe with cut followed by paste command, see manual for paste command here.

I have used this command cut -f 2,125 CADD.tsv | paste fileA but it gives me only first two columns which are 1 and 2.

cut -f 2,125 CADD.tsv | paste fileA - I think you need to put a dash at the end?

Sorry for bothering again but putting dash at the end paste only column #2, it does not give column 125

To be honest, both lines work on my machine: paste fileA <(cut -f 2,125  CADD.tsv) > myNewFile and cut -f 2,125 CADD.tsv | paste fileA - they give me 4 columns, last one called PHRED.

My first guess is that the script isn't being run by bash. What operating system are you using? What exact command do you use to launch the script? Also, both of your files have headers which are not the same so should be skipped when pasteing.

I am using Ubuntu 16.04.5. I used this command sh Cadd.rsids.sh. I do not want header just columns. When I run these commands separately, they run without error. But not in script form. :(

Your problem is that you have written a script in bash but are then running it using sh. Bash, the Bourne-again shell, is not the same as the Bourne shell (sh). What is more, on Ubuntu, sh is actually not even the Bourne shell but another minimal shell called dash. The /bin/sh on Ubuntu is a symlink to dash:
$ ls -l /bin/sh
lrwxrwxrwx 1 root root 4 Mar  1  2018 /bin/sh -> dash

Dash is a POSIX compliant shell which means it has the features defined by the POSIX specification for the sh shell. In your case, the problem is caused by the fact that the <() input redirection construct is supported by bash but not by dash. So, in order to run your script, you can do one of two things:

Make the script executable (chmod a+x cadd.rsids.sh), and then run it directly:
/path/to/cadd.rsids.sh

This is the simplest approach since you already have the shebang line (#!/bin/bash) which will cause the script to be run by bash and not sh.
Run the script with bash instead of sh:
bash /path/to/cadd.rsids.sh

Instead of
sh /path/to/cadd.rsids.sh


text:
How to paste RSIDs in CADD output
I want to paste RSIDs in CADD output as CADD does not give RSID column in its output. For this purpose I am using bedtools intersect to compare two files and have RSID column in my CADD file.
This is the command I am using.
bedtools intersect -a CADD.bed -b Input.vcf| cut -f 3

But this command is giving me error
ERROR: file CADD.bed has non positional records, which are only valid for the groupBy tool.

Is there any way to solve this error or any other way to get RSID column in CADD output file.

You might consider vcfanno.

It would be useful if you show us some line of your CADD.bed and Input.vcf and how your desired output should look like.
bedtools is complaining about non positional records in CADD.bed, which means there are lines in this file where the first three columns aren't CHROM, START, END.
I guess you will have some header lines in the beginning that aren't introduced by an #. Just delete them (or add in # to the beginning of the line) and try again.
fin swimmer