text:
CRISPR screening and systems biology
I am not sure how to tackle this one, so I’ll explain the general idea I have in mind.
Given results of a CRISPR knockout screening experiment (like in this paper) with two types of assays: control vs. treated with 5 time points per each of the assays.
Is there a way to use the screening data to perform systems biology network analysis? Something which is related to signal transduction, transcriptional regulatory, or metabolic networks?
Maybe something which uses Bayesian inference? Or something similar to differentially expressed representation of the CRISPR screening? Alternatively, a GWAS analysis similar to this paper?

do you have a list of genes that you would like to investigate through network analysis?

@gabrielet, yes I do have such a list. Even a few related such lists. But I also want to incorporate the information from the up regulation and down regulation values I got in the analysis. Not just the list of genes

What do you call screening results?
Do you have a kind of array with protein or RNA data?
If so, you can use several open-source tools available on the internet, like Enrichr or network analyst to better understand the pathways affected by your treatment.

https://www.nature.com/articles/s41467-018-06676-2

what's network analyst?

https://www.networkanalyst.ca

To investigate genes through a network based approach you can use several tools. I don't know which level of programming you have or which tool you prefer but here is a brief overview of the possibilities you have.
You can look into established databases containing information about interactions like, for instance, STRING. Or, you can build networks using the wgcna method with R.
Then, once you have your network you can analyse its properties and here you also have several tools available such as standalone, java-base software like cytoscape, R libraries such as igraph, python libraries such as networkX online tools like centiserver...
Also you can perform gene set enrichment analsysis, which means adding a further layer of biological information to you set of genes like metabolic pathways (from KEGG, for example) and molecular functions (from GeneOntology, for instance) (gsea) as here, here, and here.
Here there is also a nice tutorial about network analysis with R.


text:
TBLASTN error: Too many positional arguments
I was running tblastn using standalone blast2.9.0+ through Ubuntu Linux, to get the alignment result of a list of protein sequences for same species, against its complete genome file.
The command I was using:
~/softwares/ncbi-blast-2.9.0+/bin/tblastn -query TBLASTNstring/882.Tstring -subject ../com_gen/*882*.fasta -outfmt "7 acc qacc pident qcovs evalue sstart send sframe" -max_hsps 1 -num_alignments 1 -out 882.blastout
(Info update: The query file has more than 1 sequence but all are in fasta format; the subject file has only 1 sequence; there is only 1 file with "882" named in the folder)
And I received error from tblastn as follow:
Error:  (CArgException::eSynopsis) Too many positional arguments (1), the offending value: ../com_gen/DP_882_DeVul_Hilde.fasta

All the sequences in both my query and subject sequence files are of two lines, the header line started with ">", followed with sequences.
I have checked both query and subject sequence files using sed -n '0~2p' | grep -v -E '[[:alpha:]]' to see if there exist special character or space within the sequence line, but the output of the results return nothing.
PS: the complete genome file contained of nucleotide sequences with only ATGC (I have checked using grep -v "ATGC")
I ran the same command for other species as well and it worked fine.

Info Update: The screenshot above was example of the result for species with TaxID 196627.
I wondered what is the thing that went wrong here? And what else I can check on them?

How many files match the ../com_gen/*882*.fasta glob? Can -subject take more than one file as an argument? Also, aren't you mixing query and subject? What is TBLASTNstring/882.Tstring? That doesn't look like a fasta file, could that be the target DB name?

-subject and -query only allow one fasta file as argument

@PeterMenzel I haven't run blasts in ten years now. Do the modern versions really accept fasta files as -subject? Do they create the DB on the fly? Don't you need to first build the DB from the fasta file?

Hi @terdon I'm sorry for not being clear enough in the description. I have just updated some info in the post which is relevant regarding your question. There is only one sequence file named ../com_gen/*882*.fasta in my folder. (I was using * as the file name is quite long but certainly there is only 1 sequence file with the name of "882".) TBLASTNstring/882.Tstring is my query file which has more than one sequence but all in fasta format.

@PeterMenzel As far as I am concerned, blast can detect all sequences in a file given if they are in fasta format? Please correct me if I am wrong.

Anyway I just tried extracting one of the fasta sequences from the query file and ran it again against the genome. The same error still appear.

source: Stack Exchange

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

sort: votes


text:
How to quickly and robustly convert between mmCIF and PDB?
There is already a question on PDB/CIF to MMTF, however what is a robust way to programmatically go between PDB and CIF files?
For example I can use a python script from this gist that relies on Biopython.
id="4ckh"
wget http://www.ebi.ac.uk/pdbe/static/entry/download/$id-assembly-1.cif.gz
gunzip $id-assembly-1.cif.gz
python cif2pdb.py $id-assembly-1.cif $id.pdb

This process yields a pdb structure as shown below.

The issues are:

Requires python/Biopython to be installed.
The output isn't perfectly clean.
The gist is relatively undeveloped and hasn't got much evidence of testing.

I can't really address your issues certainly so I thought I would comment rather than answer, but I've had good results using Maxit (https://sw-tools.rcsb.org/) which is produced by the RCSB. At least for submitting to the PDB Maxit is endorsed so I would hope it's robust.

Converting a PDB file to some version of an mmCIF file is always possible provided the PDB file doesn't have any "problems". However grouping identical molecules into the same entity is non-trivial so the Biopython and BioJulia mmCIF writers (which I wrote) do not try and do this, simply assigning new entites for each molecule. These files will also be missing other things like the metadata parts of the mmCIF dictionary not involving the atoms.
In Biopython:
from Bio.PDB import PDBParser, MMCIFIO
p = PDBParser()
struc = p.get_structure("", "file.pdb")
io = MMCIFIO()
io.set_structure(struc)
io.save("file.cif")

In BioJulia:
using BioStructures
struc = read("file.pdb", PDB)
writemmcif("file.cif", struc)

Going the other way is, however, not always possible. You can run out of atom serials or chain IDs to use in the output PDB file (this is one reason for the switch to mmCIF files in general). The gist you have linked is a way to get round this.
If this is not a problem for your structure you can invert the above snippets to get a decent PDB file - see the relevant docs. Note that in both cases it will use the auth_seq_id and auth_asym_id fields to get the residue number and chain ID rather than the relabelled label versions. The test suites of Biopython and BioJulia do test these conversions pretty thoroughly.
There is also pdb-tools if you do not want to install Biopython, though both involve a simple pip install so it's not much difference.
For some speed benchmarks of various software see here and here.

hey @jgreener! Thanks for the clarification!
I'm wondering how come MMCIFIO drops metadata and such? Isn't that pretty inconvenient? let's say i need to rename (change the ids of ) certain chains in a structure and nothing else... I'll be losing all the authors' annotation, seq data etc. just by doing io.set_structure(struct), right? Are there any alternatives that you know of?

Yes, you would lose that info. The readers and writers focus on parsing the atomic structural data into a unified representation that can be modified and written back out to any file type. Developer focus has gone there rather than on writing a unified representation for the metadata. My recommendation would either be to concatenate the newly generated file onto the old metadata (beware some may be wrong now) or to use the MMCIF2Dict constructor and modify list items (writing that out preserves metadata, again some metadata may now be wrong).

that's super useful. I have played with MMCIF2Dict before and that'd actually be ideal to preserve and edit the metadata. Any tips for how to write that representation (the dict resultant from MMCIF2Dict) to a .cif file if not through MMCIFIO?

You can write a MMCIF2Dict with MMCIFIO, just use io.set_dict(d) rather than set_structure.

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


text:
How do you map nanopore fast5 files?
Is there a best practices for mapping Oxford Nanopore files to a reference?
Is there a tool that can take a tarball of fast5 files and map them directly or do they need to be converted to fastq first?

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


text:
GATK 4.1.4.0 Mutect stats output
I have adapted my pipeline to the new filtering strategies implemented in gatk 4.1.4.0 and while looking at its output I noticed that the stat file generated by mutect2 has a negative number:
callable -1.538687311E9

I did not find the documentation about this file content but in the gatk forum it is mentioned that is should represent the number of bases above a certain coverage (10 by default).
This happened for a WGS sample, while for WES I obtain a more sensible value (well, at least positive and integer, even if maybe a tad high):
callable        8.6884303E7

Has this happened to anyone else?
I compared the resulting calls obtained for the WGS with a previous version of mutect2 and they are comparable-ish (98% of the new "ones" are in common with the "old" ones), but I am still worried that something fishy is happening, expecially considering that the filtering step of the pipeline uses this stats file.

quick followup: seems to be a bug due to an integer overflow (yay Java!), GATK developers have fixed it and I'm trying it out.

This was answered by the OP -- it was a bug that has been fixed in a subsequent release, 4.1.5.0. See the release notes here: https://github.com/broadinstitute/gatk/releases/tag/4.1.5.0

Made callable sites a Long to avoid integer overflow (#6303)

Issue ticket: https://github.com/broadinstitute/gatk/issues/6302

Hey ya, this could be a useful answer if you would provide some details, like the version where the bug is fixed and a link to the changelog. Right now it provides nothing than what OP wrote.

Sure; done. I was just looking to not have it show up as unanswered anymore since the problem was effectively solved. But I get your point.


text:
Can Maven (mass spectrometry software) read mzML files?
I did not find a documentation. Can Maven (mass spectrometry software) read mzML files?

No, Maven still can read only mzXML files.


text:
How to do the reverse of bcf merge?
How can one do the reverse of bcftools merge?
To get one sample back out of the merged one.
Without holding reference sites.

bcftools has the view command which can do this. Create a text file (samples) which has the names of the samples you want, one per line, and then:
bcftools view -S samples merged.vcf.gz > new.vcf

Or, if you're running a shell that supports it, you can use input redirection for a single sample name:
bcftools view -S <(echo "sample name") merged.vcf.gz > new.vcf


text:
SRA Toolkit and lebanese data
I am trying to extract data from this:
https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR6245218
I installed SRA Toolkit, downloaded the SRR6245218 file and executed this:
fastq-dump -X 5 SRR6245218 --fasta -O ./

The result is a fasta file 900 bytes long, which does not contain the mitochondrial data that I need.
What am I doing wrong?

Your command is fine. It is working as expected. Did you intend to include -X 5 in your command? That restricts the number of spots to 5. If you drop it, you will download the entire set of 2563898 spots.

I see. I expected to get mitochondrion DNA. Not sure what I am getting here. But I already found it in fasta format somewhere else.


text:
What is the state of the art for GWAS in terms of the statistical algorithm for either Case/control and Quantitative traits?
This question was also asked on Biostars
I'm trying to understand what is the best algorithm for GWAS nowadays.
I know we have many tools available like Plink and Hail, but currently, what is the best algorithm if I won't use any them? Let's say, write down a script in R or Python from scratch. Which statistical algorithm should I use?
Is it Linear mixed models (LMMs)? I'm confused as we can have binary phenotypes (case/control) or quantitative phenotypes. LMM seems to address quantitative ones, but can it be used for case/control as well?
Actually, what is the state of the art for both/each of them? Pair-reviewed papers as references will be appreciated.
Thanks!

To clarify: you're asking what the best algorithm is, assuming you don't want to use the well-tested and validated algorithms that other people are using?

It can be whatever anyone is using. I just want to know what is the algorithm behind the scene for binary or quantitative phenotypes.

This is not a good fit for a Bioinformatics Stack Exchange question, as the question is too broad and primarily opinion-based. I'd recommend that you ask something like this on reddit. The best questions on BSE deal with specific problems that people are having. If you think your question is specific, please edit it to make that more obvious.

Sorry, but I disagree. I'm not looking for opinions, but the algorithm that, mathematically, (through statistical tests) prooves to be the state-of-the-art. I've got the answer, check it out below. Thanks anyway!

I've got the answer to my question at https://www.biostars.org/p/410466/ by chrchang522

"The main regression executed by Plink was introduced by EIGENSTRAT in ~2006; see https://www.nature.com/articles/ng1847 . This is actually straightforward to write in R/Python from scratch; the harder part is optimizing the implementation for large datasets.
The Firth regression added to Plink 2.0 to improve handling of rare variants and imbalanced binary phenotypes was motivated by https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4049324/ .
Mixed linear models provide better statistical power when you have lots of close relatives in your dataset, but are much trickier to solve; actually, this is still a significant research area. Two tools covering parts of the current state-of-the-art are SAIGE (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6119127/ ; handles imbalanced binary phenotypes, but relatively slow) and fastGWA (https://www.nature.com/articles/s41588-019-0530-8 ; great speed, but doesn't support dosage data yet and uses a misspecified model for binary phenotypes)."


text:
how to solve “dim(X) must have a positive length” at running ComBat function in R
I used ComBat() for batch effect correction in my expression data. basically, that function inputs are expression data, Batch covariate, and Model matrix for the outcome of interest and other covariates besides batch. So, I prepared all inputs based on the instruction of ComBat() function in 'sva' package. "Pheno_LMS" is my clinical data inclus batchId of my samples and "TCGA_expr_log" is my expression data. my code is below:
Pheno_LMS structure:
Sample        batchId   age  . . .
GSM71019.CEL    396     63
GSM71020.CEL    396     58
GSM71021.CEL    410     85
GSM71020.CEL    411     58
GSM71021.CEL    410     40
.
.
.
dim(Pheno_LMS)
[1] 74 37

TCGA_expr_log structure:
                   GSM71019.CEL 1      GSM71019.CEL 2        GSM71019.CEL 3  . . .
    Gene id
  ENSG00000000003     10.115170           8.628044                  8.779235
  ENSG00000000005     5.345168            5.063598                   5.113116
  ENSG00000000009     6.348024            6.663625                   6.465892
    .
    .
    .
dim(TCGA_expr_log)
[1] 42817    74

batch = Pheno_LMS$batchId
modcombat = model.matrix(~1, data=Pheno_LMS)
TCGA_expr_Co <-  ComBat(as.matrix(TCGA_expr_log),batch = batch,mod = modcombat, par.prior = TRUE,mean.only = TRUE)

but when I run that function, I got below error:
Error in apply(dat[, batch == batch_level], 1, function(x) { :
  dim(X) must have a positive length

Now, I need anybody comments about my problem. thanks
P.S: when i call batch I see:
> batch
 [1] 13 13 13 14 14 14 14 14 14 14 14 6  6  7  7  7  7  7  7  7  7  7  7  7  7  7  7  7
[29] 7  8  7  8  11 11 11 14 12 1  6  4  2  3  3  3  3  6  6  1  3  9  8  6  11 7  3  5
[57] 4  5  5  4  4  11 8  8  12 12 12 12 12 10 10 10 11 13
Levels: 1 2 3 4 5 6 7 8 9 10 11 12 13 14

dat[, batch == batch_level] is likely empty.

Thanks for @DevonRyan comment. at the end of my post, I add "batch" content. As you see, when I call batch, for each sample I see a batch value. I have 14 batches.based on your knowledge how can I solve this problem?

You need to show us something that’s reproducible. As is we need to read your mind to figure out exactly what’s going wrong.

Dear @DevonRyan, for finding problem I save "Pheno_LMS" and "TCGA_expr_log" to "my_data.RData" which you can download it from this link:
https://www.dropbox.com/s/is4b9laou305d5k/my_data.RData?dl=0

Dear @DevonRyan, I appreciate it if you check my_data.RData and guide me for solving the problem.

batch = Pheno_LMS$batchId should be batch = Pheno_LMS$BatchId, if you fix that then at least the code you showed works just fine.

Thanks @DevonRyan.In addition of your comment I found 2 batches by just 1 sample. So, I removed them from my TCGA_expr_log and now ComBat() runs perfectly.

as suggested in the comment you should verify if the command:
dat[, batch == batch_level]

returns an empty object or not. Apparently, batch and batch_level do not have anyting that match and the function tries to run over an empty dataframe.

Thanks for @gabt comment. in the end of my post, I add "batch" content. As you see, when I call batch, for each sample I see a batch value. i have 14 batches.


text:
Add new tool to galaxy
I am trying to include a new tool to the galaxy main menu, following these instructions.
tool_conf.xml doesn't exists only tool_conf.xml.main and tool_conf.xml.sample and none of them seem to be working. I make changes and nothing happens when I run galaxy.
What am I missing?
Thank you.

The .main and .sample file files usually come with default installations and may have been removed by whoever setup this Galaxy instance. You can look at the originals online and just base your own file on those: https://github.com/galaxyproject/galaxy/blob/dev/lib/galaxy/config/sample/tool_conf.xml.sample

Thanks @Devon, I removed the suffix .sample and it worked. I believe that galaxy folks could improve/update their manual.


text:
which NCBI tool is optimized to identify a species from a DNA fragment?
“The Iceman” was a man who lived 5300 years ago and whose body was recovered from an Alpine glacier in 1991.  Some fungal material was recovered from his clothing and sequenced.
Ice man : found as a frozen corpse in the Val Senales Glacier in Italy.
Thought to be carrying Fomitopsis betulina (formerly known as Piptoporus betulinus) commonly known as birch polypore, birch bracket, or razor strop) is a common Basidiomycota brown rot macrofungus growing on decaying birch wood.
Other species known : Fomes fomentarius (commonly known as the tinder fungus, false tinder fungus, hoof fungus,] tinder conk, tinder polypore or ice man fungus)
Given two DNA fragments from NCBI genbank, try to find what species these DNA fragments belongs to or what species these fragments are most closely related to.
>Z54169.1 Iceman fungal clone 'Sim C-a1' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAAACTTGGGCCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTTCCTTCTGGAGAACCTCATGCCCTTCACTGGGTGTGTTGGGGAA
CCAGGACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAAT
AATAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCG
GGGGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCAT
TTGCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAG
TCTTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGA
GAAATCAAAGTTTTTGG

>Z54155.1 Iceman fungal clone 'Sim C-a2' 18S rRNA gene (partial)
GTATATTAAAGTTGTTGCAGTTAAAAAGCTCGTAGTTGAACCTTGGGTCTGGCTGGCAGGTCCGCCTCAC
CGCGTGTACTTGTCCGGCCGGGCCTTCCTTCTGGAGAACCTCATGTCTTTTATTAGGCGTGTTGGGGACC
AGAACTTTTACTTTGAAAAAATTAGAGTGTTCAAAGCAGGCCTTTGCTCGAATACGTTAGCATGGAATAA
TAGAATAGGACGTGCGGTCCTATTTTGTTGGTTTCTAGGACCGCCGTAATGATTAATAGGGACAGTCGGG
GGCATCAGTATTCAATTGTCAGAGGTGAAATTCTTGGATTTATTGAAGACTAACTACTGCGAAAGCATTT
GCCAAGGATGTTTTCATTAATCAGTGAACGAAAGTTAGGGGATCGAAGACGATCAGATACCGTCGTAGTC
TTAACCGTAAACTATGCCGACTAGGGATCGGGCGATGTTCTTTTTCTGACTCGCTCGGCACCTTACGAGA
AATCAAAGTTTTTGG

Now I know I can just google what species of mushrooms (which is where the above background info is from). But say that I was unable to and had to just go off the DNA fragments in order to find the species.
I've tried BLASTn, but I'm not sure if it is the right tool or if I need to properly calibrate the parameters to use it.
so far, I've tried:
Under “Choose Search Set:
Database : Non-redundant protein sequences (nr)
Organisms excluded: Animalia (taxid:33208), Plantae, Prokaryota  (taxid:2), Protozoa
Organisms included: fungi  (taxid:4751)
Included: Models (XM/XP),   Uncultured/environmental sample sequences
Under Program Selection:  Optimized for Highly similar sequences (megablast)

Algorithm Parameters (General Parameters, Scoring Parameters, Filters and Masking,
PS/PHI?DELTA BLAST)

are all kept default.
```

So, what were the results? What makes you think this isn't the right approach? Did you not find any good hits?

For fungal species I think normally the ITS (internal transcribed spacer) is used for identification instead of 16s, 18s in your case.

My conclusion from the tree is the fungi appears to be genuine and not an extant contaminant (which is very easily done). It is not Fomitopsis betulina because this is division Basidiomycota. The fungi here is division Ascomycota and these are very distinct from Basidiomycota.
There have been many examples of contamination in ancient DNA work, but this doesn't appear to be one of them. The rationale is given below.

The first fragment was not interesting and described 100% identity to ascomycetes. The tree above is a simple distance tree from the second fragment of 18S which shows it is an outgroup to the ascomycetes (green terminal taxa). The ancient fungi contains 3 (I think) unique nucleotide residues that have not been sequenced in any extant ascomycetes to date and sequences are unique to any known species. The location of the 3 unique nucleotides is tightly clustered within the locus with 100% identities either side. This is significant because it demonstrates it is not a PCR mosaic artefact (yep its been done in ancient DNA work and subsequently proven artifactual). In the latter PCR mosaic artefact I'm referring to (which I wouldn't cite), was a fungi/bacterial mosaic from memory.
Given the homology of 18S in the extant ascomycetes I would say this appears to be a new genus, however you would need to be sure this was not sequencing error and this is very old DNA, which is subject to degradation even when frozen. The implication of the tree is quite unusual, it implies widespread extinction of this group of fungi in recent history (5000 years is slight in evolutionary terms) BUT it could be a sampling bias, i.e. no one has sequenced those fungi and NCBI is not representative.
I assume that you have performed phenotypic/morphological analysis on the spore to believe it is Basidiomycota.
The fragment in yellow is the query and the sequence next to it was the sequence already deposited on Genbank of this ancient isolate. You would need to look up the biology and ecology of ascomycete to understand its significance, but it is in a lineage of its own, so maybe difficult to interpret by comparison to extant species. You could concatenate the two alignments based on each fragment, but the problem with 18S is that is a hard to model and you'd need to perform secondary RNA structure modelling to assess whether these mutations would be viable.
In conclusion, it definitely looks genuine but it is not the species you think it is, in fact not even close.


text:
What software can I use to develop a population haplotype network?
I looking for software to conduct a haplotype network using cytochrome oxidase (CO1) mitochondrial sequence data to study population structure in Anopheles mosquitoes.

The software you are looking for is Network 5.0.0.3
Phylogenetic Network Software . Network © Copyright Fluxus Technology Ltd 1999–2018. Retrieved from https://www.fluxus-engineering.com/sharenet.htm
Whilst its old code it is still regularly used in pop gen for example there guys here recently used it. It will do the job, although personally I think Google Fusion tables are very cool, but that is complicated to use properly. R has some cool stuff in this context BTW, but in context I'd keep it simple (i.e. Network 5)


text:
How to convert gff to gtf?
My annotation file is in .gff format. I would like to convert it to .gtf format or to know if there is a way to directly download the annotation file in .gtf format?
I am working on sequences from the P Falicparum 3D7.

I think the question was already extensively answered here:
https://bioinformatics.stackexchange.com/questions/2068/how-to-convert-gff3-to-gtf2

Thank you for your awnser

I would suggest to use agat_convert_sp_gff2gtf.pl from AGAT because you loose information with gffread.
e.g here a gff example:
##gff-version 3
scaffold625 maker   gene    337818  343277  .   +   .   ID=CLUHARG00000005458;Name=TUBB3_2
scaffold625 maker   transcript  337818  343277  .   +   .   ID=CLUHART00000008717;Parent=CLUHARG00000005458
scaffold625 maker   CDS 337915  337971  .   +   0   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 340733  340841  .   +   0   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 341518  341628  .   +   2   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   CDS 341964  343033  .   +   2   ID=CLUHART00000008717:cds;Parent=CLUHART00000008717
scaffold625 maker   exon    337818  337971  .   +   0   ID=CLUHART00000008717:exon1;Parent=CLUHART00000008717
scaffold625 maker   exon    340733  340841  .   +   0   ID=CLUHART00000008717:exon2;Parent=CLUHART00000008717
scaffold625 maker   exon    341518  341628  .   +   2   ID=CLUHART00000008717:exon3;Parent=CLUHART00000008717
scaffold625 maker   exon    341964  343277  .   +   2   ID=CLUHART00000008717:exon4;Parent=CLUHART00000008717
scaffold625 maker   five_prime_UTR  337818  337914  .   +   .   ID=CLUHART00000008717:five_prime_utr;Parent=CLUHART00000008717
scaffold625 maker   three_prime_UTR 343034  343277  .   +   .   ID=CLUHART00000008717:three_prime_utr;Parent=CLUHART00000008717

with gffread you get:
scaffold625 maker   transcript  337818  343277  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    337818  337971  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    340733  340841  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    341518  341628  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   exon    341964  343277  .   +   .   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 337915  337971  .   +   0   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 340733  340841  .   +   0   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 341518  341628  .   +   2   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";
scaffold625 maker   CDS 341964  343033  .   +   2   transcript_id "CLUHART00000008717"; gene_id "CLUHARG00000005458";

while with AGAT you get
scaffold625 maker   gene    337818  343277  .   +   .   ID CLUHARG00000005458 ; Name TUBB3_2 ; gene_id CLUHARG00000005458
scaffold625 maker   mRNA    337818  343277  .   +   .   ID CLUHART00000008717 ; Parent CLUHARG00000005458 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    337818  337971  .   +   0   ID "CLUHART00000008717:exon1"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    340733  340841  .   +   0   ID "CLUHART00000008717:exon2"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    341518  341628  .   +   2   ID "CLUHART00000008717:exon3"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   exon    341964  343277  .   +   2   ID "CLUHART00000008717:exon4"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 337915  337971  .   +   0   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 340733  340841  .   +   0   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 341518  341628  .   +   2   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   CDS 341964  343033  .   +   2   ID "CLUHART00000008717:cds"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   five_prime_UTR  337818  337914  .   +   .   ID "CLUHART00000008717:five_prime_utr"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717
scaffold625 maker   three_prime_UTR 343034  343277  .   +   .   ID "CLUHART00000008717:three_prime_utr"  ; Parent CLUHART00000008717 ; gene_id CLUHARG00000005458 ; transcript_id CLUHART00000008717

So you can see as example that Name=TUBB3_2 from the gene feature disappeared with gffread.

You can use gffread to convert gff to gtf2, below is from the manual:

In order to see the GTF2 version of the same transcripts, the -T
  option should be added:
gffread -E annotation.gff -T -o- | more

The examples above also show that gffread can be used to convert a
  file between GTF2 and GFF3 file formats.

Thk for your awnser

Actually I started a mini-review about the tools to do such conversion.
You can find it here:
https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff_to_gtf.md
As I mentioned earlier they do not behave all the same way...

You can try bioinfokit (https://github.com/reneshbedre/bioinfokit) in python.
from bioinfokit.analys import gff

gff.gff_to_gtf(file="yourfile.gff3")
# The converted gtf file will be saved in the same working directory.
```


text:
Difference between clinical significance and statistical significance
I am working on a clinical data analysis project and am from IT domain.
While I am analyzing the data, I came to know that though some of the variables are not statistically significant, they are clinically significant.
I also read that difference in median values between case and control group for variables is used as a measure to indicate clinical significance
Can someone help me understand this with a simple example as to why difference in median between groups is considered for clinical significance?

nice and a great question

Another viewpoint could be : statistically significant implies use of sampling theory to establish presence of significant difference  between  case and control groups. The "clinical " significance seems to test for correctness or validity of measures of variables. Apparenty, the measurement theory is meant for checking validity of data - the difference  could arise on account of incorrect measurements. The non-parameteric tests -  it seems - take care of measurement errors and hence, produce perhaps more reliable results.

Thanks upvoted for the response

This is non-parametric statistics, the mean requires conformation to the normal distribution (but there are exceptions), or perhaps better put it requires a fixed relationship between mean and standard deviation. Non-parametric stats always use the median and are distribution free. Wilcoxon test/ Mann-Whitney U test  are used to replace t-test. The Kruskal-Wallis test replaces the 1-way ANOVA
Finally the Friedman test replaces 2-way ANOVA but has its critiques.  The t- and ANOVA tests require normality, F-max is one approach to determine normality.

Thanks. upvoted for the response

:-) thanks @subhash, likewise.


text:
Appropriate value to replace outliers in clinical studies
I am working on clinical study where I am focusing on binary classification problem.
There are few variables which has outliers.
I am trying to know which would be the right value to replace outliers in a clinical study?
Is it mean, median min,max or mode of a column of interest?
Currently I use min and max for values below and above 1.5*IQR threshold
I use IQR approach to identify outliers but would like to know how it is done usually in clinical studies?
Can anyone help me with this?

In my experience, it is more common to use all clinical data as-is for clinical studies. And if data is missing, either omit the sample or omit the variable with missing data.
If your classifier can't handle the wide variation commonly seen in clinical studies then you may want to use a classifier which is less impacted by outliers.

Thanks for the response. Upvoted.

Removing outliers is common practice in statistical modeling and perfectly acceptable. However, with regards 1.5 IQR I am far from certain about this approach. Normally, if you want to be conservative then 3 standard deviations (SD) denote an outlier, which is more stringent than IQR. Some use 2 SD. If the value is lower than 2 SD from the group mean it isn't an outlier, that's because 1.96 SD is the normal distribution.
Normally you would delete the outlier and it takes no further part in the calculation, assigning it to 'na' should work. The subsequent calculation will recognise an 'na' value, I think in python it is np('na') ... Np = numpy.
In this case you appear to be deleting the value and treating it as missing data. Treatment of missing data is complicated I agree 100%. Mean, median are possibilities for missing data I can't see this working for binary data. One approach is regression and in python there is a command specifically to do this (which I forget), R will have the same command and it will assign a value  for the missing data. However, you can't do this for outliers.
The other routes to remove outliers are via transformation, log transformation is the most common. Transformation is a big subject.

thanks upvoted for the response


text:
Adding new columns to this data frame
I have a dataframe like this
> dim(anno_maf)
[1] 32786   129
> head(anno_maf[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Start_Position
1:                    a       chr1       27105930
2:                    a       chr1        3342771
   End_Position
1:     27105930
2:      3342771
>

I have a list of genes like
gene name
PTEN
PTEN
MAP3K1
MALAT1
MED12
RUNX1
PIK3CA
FBXW7
ARID1A
MED12
TAB1

I want to add a column in my dataframe named COSMIC, if a gene in gene list exists in my dataframe Hugo_Symbol column, put a yes for this gene if not put a no; Finally my data frame has an extra column as COSMIC in which I have yes or no in Hugo_Symbol column based on the existence of genes in my gene list
Can you help me please?
Thank you

You can use ifelse():
anno_maf$COSMIC <- ifelse(anno_maf$Hugo_Symbol %in% your_gene_list,
                          yes = "yes",
                          no = "no")


text:
What does these terms mean?
I am working on cancer genomics failed to follow some terms even by googling
What is the difference of focal genome amplifications versus convergent amplification of given region of genome for instance fibroblast growth factors region?
Thank you for any clarification

I think these definitions still have a degree of subjectivity and I don't think they are mutually exclusive.
Convergent amplification basically means that "similar" cancers/neoplasms have "similar" amplification patterns, i.e. at around the same loci (see Myllykangas S et al.,2006 for an example).
Focal genome amplification is defined as "CNA [copy number alteration] of limited size, frequently enriched for cancer driver genes" (Krijgsmanthis et al., 2014).


text:
What are the largest fully publicly available molecular QTL datasets?
I'm familiar with the sQTL and eQTL GTEx data that can be downloaded from the
GTEx Portal, but I'm interested in gathering at least two other types of molecular QTL datasets:

Larger eQTL and/or sQTL datasets (greater overall sample size, even if performed only for blood cells)
Other functional molecular QTLs such as methylation QTLs, histone QTLs, DNA accessibility QTLs, etc.

The largest eQTL datasets I can think of is probably eQTLGen. This is a meta-analysis of many studies using >30,000 samples from whole blood. This might be less ideal if you want cell-type-specific QTL information.
For other molecular QTLs, I do not think there are large-scale datasets yet.


text:
Unable to read all alignments in stockholm format : UnicodeDecodeError: 'utf-8'
I have a file in Stockholm format and I am trying to read it using the Biopython package. The code is as follows :
import numpy as np
from Bio import AlignIO

alignments = AlignIO.parse("Rfam.seed", "stockholm")
for alignment in alignments:
    print(alignment)

It gets read for some alignments and then throws up following error :
Traceback (most recent call last):
  File "test.py", line 5, in <module>
    for alignment in alignments:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/__init__.py", line 394, in parse
    for a in i:
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/site-packages/Bio/AlignIO/StockholmIO.py", line 408, in __next__
    line = handle.readline()
  File "/Library/Frameworks/Python.framework/Versions/3.7/lib/python3.7/codecs.py", line 322, in decode
    (result, consumed) = self._buffer_decode(data, self.errors, final)
UnicodeDecodeError: 'utf-8' codec can't decode byte 0xe1 in position 2468: invalid continuation byte

How to resolve this error?

Can you not just remove the unicode character from your file? The error is due to a Unicode character (å)...

The last line indicates that this is an encoding problem. Try to see if the function you are using to read the data has an "encoding" parameter that you can set to whatever encoding your file is using. If that doesn't work, you can change the encoding of your file by:

with save as using a text editor of your choice and selecting "UTF8" as the encoding. That would be a quick work-around if you have a small number of files.

or

You can use the codecs library to check and convert file encodings within your code. Here is an example.

Warning: I haven't used Biopython in ages but encountered similar problems with Python in the past.


text:
Merging matched parts of two dataframes
I have a big data frame like this
> dim(data1)
[1] 474303    128
> head(data1[1:2,1:6])
   Tumor_Sample_Barcode Chromosome Position
1:                   A1       chr1       15772155
2:                   A1       chr1       16054813
   End_Position Reference_Allele Tumor_Seq_Allele2
1:     15772155                G                 A
2:     16054813                G                 A
>

I have another dataframe like this
Chrom   Position    Ref Alt VAF
chr1    2716062 C   A   0.323
chr1    3095112 T   G   0.365
chr1    3342771 C   G   0.123
chr1    3526376 C   T   0.365
chr1    4121415 T   G   0.225
chr1    4215706 G   A   0.125
chr1    4304014 T   G   0.286
chr1    4304016 A   C   0.302
chr1    4388809 A   C   0.135

I want to create another column named i_TumorVAF_WU in my big data frame in which VAF column comes from small data frame. Please notice that VAF values should be matched with Position in big data because small data coming from big data
Can you help please?
    > unique(merge(data1,data,x.by=data1$Position,y.by=data$Position))

Returns a weird output with losing a big part of big data frame. I just want to add corresponding VAF values to the big data by keeping all the information in big data file

Use dplyr and close some stuff to free up RAM

If you have issue with memory and dealing with large object, maybe data.table is the way to go (https://github.com/Rdatatable/data.table/wiki):
Let's take two dummy data.frame:
dfA <- data.frame(Chrom = rep("chr1",150),
                  Pos = 1:150,
                  ref = rep("A",150))
dfB <- data.frame(Chrom = rep("chr1",10),
                  Pos = 11:20,
                  Vaf = rep("VAF", 10))

The big data.frame dfA:
> head(dfA)
  Chrom Pos ref
1  chr1   1   A
2  chr1   2   A
3  chr1   3   A
4  chr1   4   A
5  chr1   5   A
6  chr1   6   A

and the small dataset:
> dfB
   Chrom Pos Vaf
1   chr1  11 VAF
2   chr1  12 VAF
3   chr1  13 VAF
4   chr1  14 VAF
5   chr1  15 VAF
6   chr1  16 VAF
7   chr1  17 VAF
8   chr1  18 VAF
9   chr1  19 VAF
10  chr1  20 VAF

Now, using data.table that is more efficient and less memory consuming for the manipulation of dataframe, you can do:
library(data.table)
setDT(dfA)
setDT(dfB)
dfA[dfB, on = c("Chrom","Pos"), vaf := i.Vaf]

And now, your data.frame dfA looks like:
> dfA[5:25]
    Chrom Pos ref  vaf
 1:  chr1   5   A <NA>
 2:  chr1   6   A <NA>
 3:  chr1   7   A <NA>
 4:  chr1   8   A <NA>
 5:  chr1   9   A <NA>
 6:  chr1  10   A <NA>
 7:  chr1  11   A  VAF
 8:  chr1  12   A  VAF
 9:  chr1  13   A  VAF
10:  chr1  14   A  VAF
11:  chr1  15   A  VAF
12:  chr1  16   A  VAF
13:  chr1  17   A  VAF
14:  chr1  18   A  VAF
15:  chr1  19   A  VAF
16:  chr1  20   A  VAF
17:  chr1  21   A <NA>
18:  chr1  22   A <NA>
19:  chr1  23   A <NA>
20:  chr1  24   A <NA>
21:  chr1  25   A <NA>
    Chrom Pos ref  vaf

As you see, we incorporate the column Vaf from the small dataframe into the big dataframe.
With your example, I will first rename the colnames of your small dataframe in order to have the perfect match with the big dataframe:
colnames(data) <- gsub("Chrom", "Chromosome", colnames(data))

then,
setDT(data1)
setDT(data)
data1[data, on = c("Chromosome","Position"), i_TumorVAF_WU := i.VAF]

Does it look what you are trying to obtain ?

Read the manuals for data.table::merge and base::merge.
This (wrong syntax):
unique(merge(data1,data, x.by = data1$Position, y.by = data$Position))

Should work fine as below:
merge(data1, data, by = "Position")

Regarding:

"losing a big part of big data frame"`

It is only returning rows when there is a match on Position, if you need to return the full dataset, then use:
merge(data1, data, by = "Position", all.x = TRUE)


text:
Changing this code in a way to work for my files
I have a bash script which extracts some information from a .vcf file but one .vcf file at each time. How I can change this script in a way to work on a bunch of .vcf files and the output is a .txt file named by the name of original .vcf
This is my script
help please

EDITED

 #!/usr/bash

#outfilename
outname=$".parsed.txt"

#Header for output file
echo -e "$outname\t$chrom\t$Pos\t$Ref\t$Alt\t$ReadCount\t$VariantAlleleCount\t$ReferenceAlleleCount\t$ReadCountControl\t$VariantAlleleCountControl\t$ReferenceAlleleCountControl\t$VAF" >> $outname ;
while read -r line ;
do

#Basic information
chrom=$(echo "$line" | sed 's/ /\t/g' | cut -f 1) #&& echo $chrom;
Pos=$(echo "$line" | sed 's/ /\t/g' |  cut -f 2) #&& echo $Pos;
Ref=$(echo "$line" | sed 's/ /\t/g' | cut -f 4)
Alt=$(echo "$line" | sed 's/ /\t/g' | cut -f 5)

#Tumor sample read, variant and reference information
ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
ReferenceAlleleCount=$(echo $ line | awk -v rc="$ReadCount" -v vac="$VariantAlleleCount" '{print rc-vac}')

#Control or Normal read, variant, reference information
ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

ReferenceAlleleCountControl=$(echo "$line" | awk -v rcc="$ReadCountControl" -v vacc="$VariantAlleleCountControl" '{print rcc-vacc}')

VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
#Print output
echo -e "$outname\t""$chrom"'\t'"$Pos"'\t'"$Ref"'\t'"$Alt"'\t'"$ReadCount"'\t'"$VariantAlleleCount"'\t'"$ReferenceAlleleCount"'\t'"$ReadCountControl"'\t'"$VariantAlleleCountControl"'\t'"$ReferenceAlleleCountControl"'\t'"$VAF" >> $outname ;

#Remove info tags from VCF
done; < <( egrep -v '#' "$1")

And the error I obtained
[fi1d18@cyan01 new]$ bash vcf_parasing.sh
vcf_parasing.sh: line 7: syntax error near unexpected token `echo'
vcf_parasing.sh: line 7: `echo -e $outname'\t' "Chrom"'\t'"Position"'\t'"Ref"'\t'"Alt"'\t'"TumorReadCount"'\t'"TumorVariantAlleleCount"'\t'"TumorReferenceAlleleCount"'\t'"NormalReadCount"'\t'"NormalVariantAlleleCount"'\t'"NormalReferenceAlleleCount"'\t'"VAF" > $outname'
[fi1d18@cyan01 new]$

My script
https://www.dropbox.com/s/di49mij2a43i43s/vcf_parasing.sh?dl=0
Two samples of .vcf
https://www.dropbox.com/s/ryk7q0aybeo5cy2/LP6008202-DNA_B03_vs_LP6008201-DNA_B03.snp.pass.vcf?dl=0
https://www.dropbox.com/s/u44s6d54r5d57bu/LP6008336-DNA_H02_vs_LP6008333-DNA_H02.snp.pass.vcf?dl=0

Doing this sort of thing with a bash loop over every line is a really, really bad idea. It will be slower than just about any other method, the syntax is hard to read (and I say this as someone with considerable experience with shell scripts) and the whole thing will be more painful than needed. If you could [edit] your question and explain exactly what information you are trying to extract, I'm sure we could help you find a better way.

Thank you @terdon. This script extract some fields from .vcf file output of Strelka mutation caller. the problem is that except for the some columns, the location of the other information is of easily detecting

Note that you don't need all of those quotes for the echo. Just use echo -e "$outname\t$chrom\t$Pos\t$Ref\t$Alt\t$ReadCount\t$VariantAlleleCount\t$ReferenceAlleleCount\t$ReadCountControl\t$VariantAlleleCountControl\t$ReferenceAlleleCountControl\t$VAF" >> $outname ;

Using shell loops for text practice is considered bad practice. It is exceedingly slow, the syntax is complicated so it's very easy to get it wrong and it's just painful. The shell isn't designed as a proper scripting language, so while it can be (ab)used that way, you really should avoid it whenever possible.
Here's a Perl script that should do exactly what your bash script does but will be several times faster and can deal with multiple input files:
#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile) or
      die("Couldn't open '$outFile' for writing: $!");
  open(my $inFileHandle, '<', $inFile) or
      die("Couldn't open '$inFile' for writing: $!");

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    my @infoFields = split(/;/, $info);
    my $readCount = $infoFields[12];
    $readCount =~ s/ReadCount=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCount = $infoFields[25];
    $variantAlleleCount =~ s/VariantAlleleCount=//;
    #/#Useless comment to fix SE syntax highlighting
    my $referenceAlleleCount= $readCount - $variantAlleleCount;

    my $readCountControl = $infoFields[13];
    $readCountControl =~ s/ReadCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $variantAlleleCountControl = $infoFields[26];
    $variantAlleleCountControl =~ s/VariantAlleleCountControl=//;
    #/#Useless comment to fix SE syntax highlighting

    my $referenceAlleleCountControl = $readCountControl - $variantAlleleCountControl;
    my $VAF = $infoFields[27];
    $VAF =~ s/VariantAlleleFrequency=//;
    #/#Useless comment to fix SE syntax highlighting

    my $outString = join("\t", $outFileName, $chrom, $pos, $ref, $alt, $readCount,
                         $variantAlleleCount, $referenceAlleleCount, $readCountControl,
                         $variantAlleleCountControl, $referenceAlleleCountControl, $VAF);
    print $outFileHandle "$outString\n";
  }
  close($outFileHandle);
  close($inFileHandle);
}

Save that file as parser.pl and then run it on all your vcf files with:
perl parser.pl /path/to/vcfs/*vcf

That should create a .parsed.txt file for each input file in the target directory.

If you just have to use the shell for this, you can try with this improved version of your script:
#!/usr/bin/env bash

for file in "$@"; do
  outFile="${file}.parsed.txt"
  outFileName=$(basename "$outFile")
  (
    printf '%s\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t' "$outFileName"
    printf 'TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount\t'
    printf 'NormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n'
  ) > "$outFile"

  echo "Reading $file"
  grep -v '^#' "$file" |
    while read -r line; do
      read chrom Pos Ref Alt rest < <(awk -F"\t" '{print $1,$2,$4,$5}' <<<"$line")
      ReadCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 13 | sed 's/ReadCount=//' )
      VariantAlleleCount=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 26| sed 's/VariantAlleleCount=//')
      ReferenceAlleleCount=$(($ReadCount - $VariantAlleleCount))

      ReadCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 14 | sed 's/ReadCountControl=//')
      VariantAlleleCountControl=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 27 | sed 's/VariantAlleleCountControl=//')

      ReferenceAlleleCountControl=$(($ReadCountControl - $VariantAlleleCountControl))

      VAF=$(echo "$line" | cut -f 8 | sed 's/;/\t/g' | cut -f 28 | sed 's/VariantAlleleFrequency=//')
      printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$outFileName" "$chrom" "$Pos" "$Ref" "$Alt" \
             "$ReadCount" "$VariantAlleleCount" "$ReferenceAlleleCount" "$ReadCountControl" "$VariantAlleleCountControl" \
             "$ReferenceAlleleCountControl" "$VAF" >> "$outFile"
    done
done

Just to illustrate what I mean about the shell being slow, I ran both scripts above on the example files you provided:
$ time foo.pl *vcf 2>/dev/null

real    0m0.051s
user    0m0.047s
sys     0m0.003s

$ time bar.sh *vcf 2>/dev/null

real    2m28.842s
user    3m33.648s
sys     1m12.456s

As you can see above, the perl script took less than a second, while the bash script took almost two and a half minutes!

Finally, if you find the perl script too long, here's the same thing in a more condensed style:
#!/usr/bin/env perl

## Iterate over all files given at the command line
foreach my $inFile (@ARGV) {
  print STDERR "Processing $inFile...\n";
  my $outFile = "$inFile.parsed.txt";
  $outFile =~ m|([^/]+)$|;
  #/#Useless comment to fix SE syntax highlighting
  my $outFileName = $1;
  open(my $outFileHandle, '>', $outFile);
  open(my $inFileHandle, '<', $inFile);

  ## Print header
  print $outFileHandle "$outFileName\tChrom\tPosition\tRef\tAlt\tTumorReadCount\t" .
      "TumorVariantAlleleCount\tTumorReferenceAlleleCount\tNormalReadCount" .
          "\tNormalVariantAlleleCount\tNormalReferenceAlleleCount\tVAF\n";
  ## Read all lines of the input file
  while (<$inFileHandle>) {
    ## Skip headers
    next if /^#/;
    ## Get VCF fields
    my ($chrom, $pos, $name, $ref, $alt, $qual,
        $filter, $info, $format, $values) = split(/\t/);
    $info=~s/[^;]+=//g;
    my @infoFields = split(/;/, $info);
    print $outFileHandle (join "\t", $outFileName, $chrom, $pos, $ref, $alt,
                               $infoFields[12],$infoFields[25],
                               $infoFields[12]-$infoFields[25],
                               $infoFields[13],$infoFields[26],
                               $infoFields[13]-$infoFields[26],
                               $infoFields[27]) . "\n";
  }
}

Oh my God, this is amazing. I am wondering when some people put their time and efforts to helping the other people   this nicely . The perl done that in a second

Sorry @terdon your code works on Single Nucleotide variants (SNVs) calls very well but on INDELs returns empty fields that may be because for example the position of ReferenceAlleleCount in SNV and INDEL are not the same. This is link to one of my INDEL files coming from the same Strelka variant caller https://sotonac-my.sharepoint.com/personal/fi1d18_soton_ac_uk/Documents/Microsoft%20Teams%20Chat%20Files/LP6008202-DNA_B03_vs_LP6008201-DNA_B03.passed.somatic.indels.vcf I am not sure even I can post this comment or not because surely I see myself too annoying

@Angel I can't access that file. It's probably best if you ask a new question, and post some lines from the bad files. Or, if the position is different, just try modifying the script and changing the numbers in the various $infofields[N].

Thank you @terdon, I have asked this as a new question here https://bioinformatics.stackexchange.com/questions/11074/looping-over-several-files-in-bash  as I found the script specific to INDELs differs to what works for SNVs. If possible have a look please; thank you in advance

@Angel you will need to explain what information you need. Since your script doesn't work, that won't help us understand your objective. So, show us your expected output from a few lines of example input.

Thank you so much for considering my issue. Actually all I want is Variant Allele Frequency (VAF) which is obtained by Depth of Alternative Allele / Total Depth or AD/DP. I have addressed this in my asked main post.

Here is a quick and kind of dirty solution with minimum intervention to the code:

Add $outname'\t' right after echo -e to two two lines, the lines after #Header for output file and #Print output. The former will alter your header and the latter your actual output: you will have a first column with file name information.
Then you can run this script in a for loop to generate multiple output files which then can be merged into a single file with cat.

Note that with this workaround, you will have multiple headers in the merged file. If you need a clean output, you will need to put your while loop in a for loop.

Thank you so much @haci . I edited the script as you suggested adding for loop to that in my main post but I obtained the same error posted in my main post

What happens when you get rid of the whitespace in the $outname'\t' $chrom bit (remember there are two of these for the two lines that were altered from the original code) and the quotes as @terdon suggested in their comment?

Sorry I tried all of your comments but I failed I have attached my original script and two vcf files here maybe it is better than pasting the code here


text:
How can I not show insertions in the Integrative Genome Viewer (IGV)?
I am using IGV 2.5.2 and would not like to see the insertions  in my aligned reads. How can I do that? I have managed to remove mismatched bases because there is an option to do so when I right-click on the window ("Show mismatched bases"). Now I would like to do the same with these purple insertions.

3rd January 2020 update

Thanks to a member's response, I found some parameter in the preferences that might be helpful to hide the indels. However, I have tried different numbers from 0 to 1000 and the indels are still shown. What numbers should I put to the parameters "Show indel threshold (bases)" and should I click the box or not?

Check the box and try 1 as the threshold.

Preferences -> Alignments -> Hide indels < 'x' bases.

Credit to Pierre Lindenbaum


text:
Are there any "canonical" sources of C. elegans connectome?
I found that there are a couple of projects that allow the user to simulate the neural system of C elegans, for example OpenWorm one update being here.
I don't know how to navigate sources in bioinformatics and thus struggle to evaluate their quality.
Question Are there any fundamental / most widely accepted data sources for raw data of C Elegans connectome or electron microscope scans of this organism, similar to Human Genome Project? If not, what are the ways to derive the data?

Answer from @terdon, converted from comment (please improve if possible):
Have you checked on WormBase? That's the main repository for c. elegans data. I don't know that they have connectome data, but if they exist, that's where they will be stored.
[see here for more information on the purposes of comments in the StackExchange network]


text:
Looping over several files in bash
I have a bash script:
I am wondering how I can change this script to loop over a bunch of .vcf files and give output .txt with the name of corresponding .vcf
I tried changes done in similar script in  Changing this code in a way to work for my files but I just  messed everything up

If we save the script pasted in the main post as a sh file and we have some .vcf files in a folder, by this line we can iterate over vcf files to extract what mentioned in the script returning .txt output by the name of corresponding vcf
for file in *.vcf ; do bash indel_vcf_parasing.sh $file ; done

Can you provide some context to future readers about this answer

This is generally the method I use for scaling. Write a program to do 'one', then write wrappers, or use one-liners (as this answer shows) to do many.
You will need to be sure that successive runs do not step on each other (if using a single directory for all runs), or make a new directory for each run if necessary.

Please edit your answer to include the extra details so the answer is complete and a user doesn't have to read all the comments to understand.

Error means that quotations are missing. I am not going to review this chunk but at first glance check
printf 'TumorAltAlleleCount\tTumorAltAlleleFrequency\t'\n'

There are three quotation marks so one is not matched. I suggest you develop code within a good IDE that will help you debug. Please google around for potential solutions. I personally do R and Shell scripting in RStudio. It differently colors quoted and unquoted code so you immediately see if you missed a closing quotation somewhere.


text:
Infer the new (Ban et. al) ribosomal nomenlature (ex. uL53 ) from the ribosomal protein's sequence

I'm trying to programmatically construct a name for each protein subchain in any ribosome from Uniprot in accordance with Ban et. al's 2014 proposal (excerpt given below) using PDB's and Uniprot's APIs.
This is to be a part of a bigger ribosomal analysis suite released next year and we would like to standardize the names as much as possible.
Ban et. al's proposal: [https://bangroup.ethz.ch/research/nomenclature-of-ribosomal-proteins.html]

I'm however not remotely a biologist (I work in compilers) and would really appreciate some clarification on the way nomenclature works across domains of life and how homology plays into it.
This is my naive approach thus far:

I get a ribosome from PDB let's say (3J9M or 5MYJ or any other for that matter) and split it into protein subchains. So on the order of 20-80 proteins.
For each protein subchain, I can get all the information in the world on it from Uniprot: names, sequence, organisms containing it, publications, synonymous names, you name it...
Right now (again, naively) for those proteins that don't already have a new name assigned to them, I take name-synonyms that PDB has plenty of (from various publications, I assume) and scan them with regexes (ex./[LS]\d{1,2}/g) for things like "L8" or "S15".
Thus, for each subchain, I end up with something like "L15" or "S28" in multiple copies (depending on the abundance of synonymous entries in PDB).
The rut i'm in right now is this: having something like "L15" to match it to a new "Ban-nomenclature" name: they provide the nomenclature look-up tables in terms of these three "historical" naming conventions(i.e "human", "yeast", "bacteria") and a taxonomic range while all the I have is the ribosomal structure(ex `` from which I got the subchain. I'm sure there's plenty of information that I can get

Could somebody elucidate to me which assumptions I can and cannot make in this field and how exactly to use the taxonomic range?
Say, (i) can I assume that all archean names are the same as bacterial homologs?
If I draw a eukaryotic cytoplasmic ribosome, do I match its subchains against yeast column or human column (since both are eukaryotic)?

Somebody suggested that I start grouping these subchains based on sequence, not on a domain and I do have access to each protein's sequence, but I am not sure how to go about it. I was thinking to get each of Ban's new names' UniRef90 cluster, let's say, and then start checking which cluster a given subchain belongs to.
But I'm obviously guessing at this point. Any advice or pointers are very appreciated.

Uniprot90 is too limited, you'll likely get better results from  pfam id (domain family annotation) —that is, part of your proteins will be annotated in a group and fusion protein (if any) will have two annotations. Alternatively, (original) COG groupings are really handy for this but aren't available from Uniprot. Then there's always the option of making your own psi-Blast motif on the NCBI website (there's a tactic to it for best results though) and using local blast with it.

@MatteoFerla thanks a ton, this is super informative! A bit too dense, actually.
Rolling with your first suggestion: as i understand it, a "new name" (ex. "uS15") would be represented by a Pfam-group of homologs and then for a given protein with an old-name i would track it to its Pfam group?

Yes. Correct. What you are trying to do is find what genes are in what cluster, which is straightforward. You have given a researched questions as per guidelines, however, there is honestly too much to read for my lazy self (I skimmed and do not care to read the paper, which from the title reminds me of the XKCD about standard), so I am bound to be missing something. But here are my thoughts of 4 pitfalls: (cont'd)

(A) you may have two protein clusters that have a single PFAM, (B) you have gene fusions (one protein, 2 domains of interest), (C) the a names are meant to cover analogues, i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural), (D) you have a very small protein which is poorly conserved, (E) a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).

This is a copy-paste nearly verbatim of comments so that the question has an answer
What you are trying to do is find what genes are in what homology cluster. This is common problem and there are many solution each with some issues.
Uniprot90 is indeed a cluster of homologues but it is too limited. Whereas you require clusters that span all the domains (universal).
COG
The nearly original grouping, is COG (cluster of orthologous genes). It is really handy for this but aren't really available from Uniprot and NCBI (its curator) has been phasing COGs out. In Uniprot records, they are patchily present as its derivative eggNOG (which has varying degrees of "scope", LUCA being the universal one):
<dbReference type="eggNOG" id="COG0451">
<property type="taxonomic scope" value="LUCA"/>
</dbReference>

Pfam
Another option that is well mapped and always present when possible in Uniprot entries is Pfam ids. These are domain family annotations and are also good option. It is a domain fold group. So domains within protein with homologous folds are grouped into one. If a protein has two known domains, it will have two Pfam entries —repeat proteins have loads, but most protein have one domain. In Uniprot XML you will find these are:
<dbReference type="Pfam" id="PF01370">

Similarly there is InterPro, which is more narrow than Pfam, but generally overlaps.
DIY
Then there's always the option of making your own psi-Blast motif on the NCBI website (setting high number of hits, setting the database to refseq or even PDB and doing a few iterations), saving the motif pattern (PSSM) and using local blast with it.
Caveats
In general you may however encounter these pitfalls:

you may have two protein clusters that share a single id
you have gene fusions (one protein, 2 domains of interest). In Pfam the protein will have two ids, with COGs or similar, you'd have only an anamalous length to go by.
you may have wished to find analogues (i.e. two protein that arose independently that have the same role (a lot of ribosomal protein are simply structural). No homology clustering scheme can find these.  Say a bacterium has evolved to have a different protein with that role, yet as it has never been studied it is a domain of unknown function (guilt by association is the common strategy, but that is beyond your scope here).
you have a very small protein which is poorly conserved


text:
Low Fraction of usable antibody reads in CiteSeq
we performed a combined gene expression and CiteSeq experiment with the 10x VDJ kit and 20 conjugated antibodies and sequenced on hiseq. I used cellranger to process the sequencing output. The cellranger summary shows overall good values except for a low percentage of "Antibody Reads Usable" and a low fraction of "Antibody Reads in Cells". How is it possible that so few antibody reads are usable? Does this mean that most antibodies were not bound to the cell surface?

We have also seen a very similar issue and are also investigating why. Possibly of note, I tried forcing the cellranger program to assign a very large number of cells (10x the actual number of cells we had) using the --force cells option. This appeared to increase the % of usable Ab reads and Ab reads per cell metrics, suggesting that these reads may be associated with barcodes that have relatively few UMIs?

Same issue here!
I do not understand how could the fraction of usable antibody reads be so low if the valid ADT UMIs, cell barcodes and antibody barcodes are all close to 90%


text:
Moving file based on their names
I have a list of vcf files; I also have a list of names in a txt file like
LP6005409-DNA_F01
LP2000325-DNA_A01
LP6005409-DNA_E02
LP6005500-DNA_C03

What I have in this txt files are part of name of each vcf files
I want to move .vcf files having these in their name to a separate folder
I done
[fi1d18@cyan01 TRG]$ for i in *.vcf
> do if grep -q $i 1.txt;
> cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
> done
-bash: syntax error near unexpected token `done'
[fi1d18@cyan01 TRG]$

How I can complete this task?
Thank you

First of all you are getting an error because you are missing a then in the third line of your code and then a fi to close the if. It should be:
for i in *.vcf
do if grep -q $i 1.txt; then
    cp *$out* /temp/hgig/fi1d18/TRG45/snp/snp/TRG/pre/
fi
done

And I guess the one-liner below achieves what you are trying to do:
while IFS= read -r line; do mv "$line"*.vcf some-directory; done < file-list.txt

for i in $(cat text.txt); do cp -v *${i}*vcf mydir; done


text:
mapping and de novo assembly of plant without reference genome
I'm studying on a plant that has no reference genome and only has one scaffold assembly and one gff3 annotation file. Can I create an index with the same assembly and gff3 in STAR and do the mapping? If I want to do de novo assembly using the Galaxy Can I import the links of 20 fastq files from the RNA-seq of the samples together to construct the Reference Genome with Trinity? If this is possible, please help me import the links of these files together into Galaxy and make the reference genome.
Thanks a lot

You have multiple questions, I will try to answer one by one:

On creating an index with STAR: As long as you provide a fasta file and its corresponding gff, you will be able to generate an index.
On creating a reference with Trinity: You have mentioned "to construct the Reference Genome with Trinity", however, with your RNA-seq data and Trinity, you can construct a de novo transcriptome and not a genome.
Uploading files to Galaxy: I suggest to do so via FTP using FileZilla or something similar. The transfer can take some time if your files are large though.

If I understand your questions correctly, you are performing RNA-seq and you are concerned about the completeness of the genome of your model organism. In that case, I would advise to bring together a de novo transcriptome from your RNA-seq data (with Trinity), use one of the quasi-mapping methods (with a tool like Salmon or kallisto).


text:
bedtool intersect to compare .vcf file to reference genome
I have a .vcf file generated using samtools. I want to know mutations in specific genes by comparing .vcf files with a reference genome. I cannot use the Ensembl VEP tool as my reference genome is not listed there. Can you please help?

Using SnpEff (http://snpeff.sourceforge.net/SnpEff_manual.html#databases) you can build custom annotation databases using the gff or gtf file. VEP also offers similar functionality https://asia.ensembl.org/info/docs/tools/vep/script/vep_custom.html.

If you want to filter a VCF file, you might want to give VCFtools a try.
Here is the documentation : VCFtools Manual
I am not a specialist of this tool but you could do as the following. First build a bedfile with the chromosome names and positions of your genes (start and stop) such as this :
chr   start   stop   (name)   (score)   (strand)

Then it should work with this command line :
vcftools --vcf file.vcf --bed positions.bed --out filtered.vcf

The output will be a VCF filtered on the provided regions

thanks, I will try that!


text:
Single-cell sequencing dataset has too many barcodes
I am analyzing a single-cell sequencing dataset from the website 10xgenomics, with 2000 cells. It is a BAM file and I am trying to obtain the individual cells per sample. I used the command
samtools view mtdnaAsorted.bam | grep CB:Z: | sed 's/.CB:Z:([ACGT]).*/\1/' | sort | uniq -c > reads_per_barcode_mtdnaA

to get a list of the individual bar codes and their reads, but there are 100,000 unique barcodes even with only 2000 cells. Does anyone know what could be the problem?
I know that there could be noise, but this is a problem for me because I'm trying to analyze coverage per cell, and if there are a lot of cells that contribute towards the total number of reads, the net coverage per cell is going to be really low

Cross-posted https://www.biostars.org/p/415116/

The BAM file you downloaded has all of the detected 10X barcodes instead of those that represent cells.  In a given 10X experiment the number of input barcodes vastly outnumbers the input cell count.  After the reaction some of barcoded beads are invariably sequenced because they are encapsulated with ambient RNA, dead cells etc. After processing the BAM file Cell Ranger infers which cells likely represent cells and generates gene expression matrices from these files.
You would need to download one of the output files that contains the list of passed barcodes such as the Gene / cell matrix (filtered) and then use the list of passed barocode sequences to filter your BAM file. Alternatively you could approximately select the same cells by selecting the top 2,000 barcodes with the highest numbers of reads.


text:
Read-informed statistical phasing
I'm wondering if there are any tools out there that perform statistical phasing of genetic data (e.g. Eagle) that take read-based phasing from GATK into account.  I've found this paper which sounds like what I am looking for, but that tool was developed specifically for Hi-C data, and I only have exomes.
Thanks

shapeit4 is currently (probably) the most accurate and quickest phasing algorithm which can account for sequencing reads in the model. They suggest you extract the phase information from the reads using WhatsHap first however. The guide for WhatsHap can be found here and has some info on using data processed using GATK.


text:
Genetic Genealogy related to SNP variants and INDEL variants
Background:
I am looking for some advice on the use of SNP variants versus INDEL variants for determining Genetic Genealogy, specifically for determining accurately the "START" and "END" of shared DNA segments (shared between two DNA testers' results). The test data used is RAW DNA from direct to consumer DNA testing companies such as AncestryDNA and 23andMe.
Question:
For the purpose of Genetic Genealogy, are the INDEL variants as useful as the SNP variants? Should the INDEL variants be ignored? Or, alternatively, do the INDEL variants hold a greater significance than the SNP variants?
I hope that using INDEL variants would increase the data in the comparison algorithm, thus increasing the certainty and resolution when determining each shared segment.

InDels are just as informative, but are less reliable to call. Errors should be averaged out in this kind of analysis though (so long as they don't predominate).
Both SNPs and InDels are valid data points to consider when comparing triads, but identifying IBD segments will only ever be probabilistic.

first of all there's no way to accurately determine the start/end position of shared DNA segments when using microarray data from DTC's like 23andMe and Ancestry.
INDEL's are indeed important when it comes to identify shared segments though I'm not sure how often the change over the course of several generations.
My mum's sister has an INDEL that both my mum and her brother didn't inherit, it will be inherited by her children (my 1st cousins).
Unfortunately the microarray from 23andMe doesn't give you any detailed information of than DI, II and II but not how long the deletion or insert is.
To determine such it requires to use WGS data.

Welcome, for the answer to be valid you would need to share the putative indel location. The familial relationships could/should be anonymised but the biological information can't.

Very informative. Thank you @a4xrbj1. I found that the length of the deletion can be read from the [link]https://www.ncbi.nlm.nih.gov/snp/ webpages, for example: [link]https://www.ncbi.nlm.nih.gov/snp/rs267607236.

@MichaelG. - I don't get your comment. This isn't a question about a specific INDEL, therefore it's a general answer. No idea why your asking for my aunt's specific INDEL at one chromosome and how that should help the OP.


text:
Why don't all reads have adaptors?
I got NGS reads back from sequencing platform. I check for adaptors and trimmed them. But I realized only a fraction (eg 30%) have adaptors... why not all of them?
Thanks

Things like this might depend on your specific library prep, but in general: sequencing starts at the end of the adapter, not before it. You will only see adapters if you sequence through the entire fragment into the adapters on the other side. If your fragment is long enough you won't see any adapter, and of course, that's the desirable outcome of your library prep.


text:
Any user friendly way to find rare mutations in whole genome raw?
Is there any user friendly way to find rare mutations in the individual human whole genome sequencing raw data? (from Dante, 30x coverage).
To be more specific, I want to find mutations from this paper:
https://docs.google.com/document/d/1EkRMuD6J0-zyMY3MegKhz7Th0hzhzVXg6DWg2Hm6l2Y/edit?usp=sharing
(it's very short, less than one page).
But I'm confused. In their paper they just talk about genes (RAD21, B3GAT2, SMC3, SCN11A , SCN5A, SCN9A, SCN10A, SCN11A, TRPA1), but not mutations? Can we find that diseased-mutations which they talking about in genes, which their list? Or is there not enough data in the paper/study for this?
And if there is not enough data, which data I need to request from authors?
Or those genes have kind a "gold standard sequences" and if that sequence different from standard - there is "diseased" gene?
(I have only very basic genetics and bioinformatics knowledges)

Cross-posted: https://www.biostars.org/p/415417/

@ATpoint do you have a script set up to systematically find when users cross-post to BioStar or Reddit? (Just wondering how common this trend is)

No, there is no system to it. Purely by chance that a title that looks familiar catches the eye. I would not say it is a trend and only happens occasionally. This toplevel question is also at https://www.reddit.com/r/bioinformatics/comments/elivcf/any_user_friendly_way_to_find_rare_mutations_in/

That isn't a paper, it's just an abstract. You can read the actual paper here: https://doi.org/10.1016/S0016-5085(17)30765-5 but it isn't available for free. Presumably, the paper will give the actual variants found.

You need to run a SNP finding pipeline, filter for SNPs in those genes and then run a SNP annotation pipeline. (You could annotate before you filter but it will take longer).
SNP finding usually has three steps, 1) read alignment to reference, 2) SNP calling and 3) SNP filtering (by quality).
You then bring in gene annotation data to filter SNPs to just those genes.
You then run SNP annotation using tools like VEP, SIFT, PolyPhen, etc.
People use tools like GATK to do this.
Hopefully that's enough information to get you started.


text:
include a glimma interface in a shiny app
I am trying to code a shiny app for RNA-Seq data analysis.
I would like to include glimma interactive plots in it.
However, in my current interface, clicking the action button Glimma plot opens the glimma interface in a new page of my web browser and not in the shiny UI.
I use shiny, glimma and EdgeR. Here is my code:
ui <- fluidPage(
  actionButton(inputId="run.glimma", label="Glimma plot"),
  htmlOutput("glimma")

server <- function(input, output, session){
  observeEvent(input$run.glimma, {
    glMDSPlot(dgeObj.norm(), labels=rownames(dgeObj.norm()$samples), groups=dgeObj.norm()$samples[,as.character(input$mds.grouping.feature)])
  }),

  output$glimma <- renderUI({
      includeHTML("glimma-plots/MDS-Plot.html")
    })
}

I only indicated the part of the code that I use to generate the glimma MDS plot.
Is there a way to display the glimma plots in the shiny app UI and not in a new webpage ?
Thank you

I have come across this page today with same query. After visiting here and there, finally following approach helped to resolve this issue. I must confess that I am just a beginner to R Shiny.
Lets say we have ui.R, server.R and global.R.
Inside global.R, add following code:
glimmaMDSOutput <- function(outputId, width = "100%", height = "400px") {
        htmlwidgets::shinyWidgetOutput(outputId, 'glimmaMDS', width, height, package = 'Glimma')
}

renderglimmaMDS <- function(expr, env = parent.frame(), quoted = FALSE) {
        if (!quoted) { expr <- substitute(expr) }
        htmlwidgets::shinyRenderWidget(expr, glimmaMDSOutput, env, quoted = TRUE)
}

Inside ui.R, add following:
glimmaMDSOutput(outputId = 'MDSplot')

**Inside server.R, add following:**
output$MDSplot <- renderglimmaMDS({
                Glimma::glimmaMDS(x = .., labels = .., groups = .., top = .., gene.selection = "pairwise", launch = FALSE, html = NULL, width = 1000, height = 500);
})


text:
Why I am getting different gene lengths for the same genes in different samples with rsem-calculate-expression?
I am running the STAR-RSEM pipeline, I have a simple shell script(*), which I use to loop over fasta files corresponding to different samples. From the output, I am interested in the *.genes.results files, for which the columns are:
> names(my_df1)
[1] "gene_id" "transcript_id(s)" "length" "effective_length" "expected_count" "TPM" "FPKM"

To my surprise, whereas the gene_id and transcript_id(s) columns are the same, I see different length and effective_length values for the same genes in different sample results (omitting the transcript_id(s) column for formatting purposes):
> my_df1 <- fread(paste0("alignment/", my_files[1]))
> head(my_df1[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count   TPM
1: ENSG00000000003.14 2292.23          2242.23            978 26.69
2:  ENSG00000000005.5 1339.00          1289.00             11  0.52
3: ENSG00000000419.12 1024.45           974.45            469 29.45
4: ENSG00000000457.13 4204.00          4154.00            241  3.55
5: ENSG00000000460.16 2953.40          2903.40             62  1.31
>
> my_df2 <- fread(paste0("alignment/", my_files[2]))
> head(my_df2[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count    TPM
1: ENSG00000000003.14 2260.38          2210.38       10052.00 117.78
2:  ENSG00000000005.5 1269.22          1219.22         132.00   2.80
3: ENSG00000000419.12 1058.13          1008.13        4386.00 112.68
4: ENSG00000000457.13 4175.05          4125.05         873.00   5.48
5: ENSG00000000460.16 2783.92          2733.92         780.01   7.39
>
> my_df3 <- fread(paste0("alignment/", my_files[3]))
> head(my_df3[1:5,c(1,3:6)])
              gene_id  length effective_length expected_count  TPM
1: ENSG00000000003.14 2662.59          2612.59             49 2.90
2:  ENSG00000000005.5  940.50           890.50              0 0.00
3: ENSG00000000419.12  956.87           906.87             58 9.89
4: ENSG00000000457.13 6364.00          6314.00             66 1.62
5: ENSG00000000460.16 2661.00          2611.00             15 0.89

The length and effective_length are defined per gene (and its transcripts) and should not differ across samples. Below is from the manual:

A gene's 'length' and 'effective_length' are defined as the weighted
  average of its transcripts' lengths and effective lengths (weighted by
  'IsoPct')

I would appreciate any hints pointing out what I am missing.
* The script I am using:
for file in ...; \
do rsem-calculate-expression \
-p 8 \
--star \
# --output-genome-bam \
--star-path path-to-star-executable \
--star-gzipped-read-file \
--forward-prob=0 \
$file \
path-to-genome \ # should end with the string used for the "reference_name" parameter of "rsem-prepare-reference"
${file:11:7}; \ # sample prefix
done

EDIT: For those who will be interested, this entry from the RSEM Users Google Group demonstrates how to calculate gene length from weighted transcript length.

If there are differences in the expression of individual transcripts between samples then both the length and effective length will differ between samples. It is commonly the case that the ratio of longer and short transcript expression varies a bit between samples, which gives rise to this.


text:
What is the meaning of the "*" character in bwa fastmap's output?
I am mapping kmers back to a few bacterial genomes using bwa fastmap:
bwa fastmap -l 9 ref.fasta kmers.fasta > out.fastmap
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[main] Version: 0.7.17-r1188
[main] CMD: bwa fastmap -l 9 ref.fasta kmers.fasta
[main] Real time: 0.168 sec; CPU: 0.157 sec

The output looks like this:
SQ      AAAAGTAGTAAGCAGGAAGACAACACGGTTG 31
EM      0       11      1       contig1:+2041368
EM      3       13      2       contig1:+1875252        contig1:-2474744
EM      4       15      1       contig1:+3779779
EM      5       16      1       contig1:-3253348
EM      6       17      2       contig1:+938710 contig1:+1066682
EM      7       20      1       contig1:-4912797
EM      10      21      4       contig1:+4803473        contig1:-2830252        contig1:+4495150        contig1:-4907283
EM      11      23      1       contig1:-3148132
EM      13      24      1       contig1:-4651172
EM      14      25      3       contig1:-2142223        contig1:+4994474        contig1:+873066
EM      16      26      4       contig1:+156775 contig1:+27749  contig1:+2207492        contig1:-1340811
EM      17      29      1       contig1:+3523989
EM      19      30      24      *

EM      20      31      3       contig1:-2533354        contig1:-208660 contig1:-1080177
//

My question is: what is the meaning of the * char in the row before last? Why is bwa fastmap not reporting all the hits as in the other lines? I could not find a page that explains this output anywhere.

Just starting to work with bwa fastmap and indeed the documentation is very sparse...
This is what we get from running the program without arguments:
bwa fastmap

Usage:   bwa fastmap [options] <idxbase> <in.fq>

Options: -l INT    min SMEM length to output [17]
         -w INT    max interval size to find coordiantes [20]
         -i INT    min SMEM interval size [1]
         -L INT    max MEM length [2147483647]
         -I INT    stop if MEM is longer than -l with a size less than INT [0]

Which is a bit cryptic as is...
The output does not really follow SAM specification and I'm pretty sure the '*' means does not mean unmapped (in fact an unmapped Exact Match doesn't make any sense does it?)
In fact, the '*' means that there are too many matches to print the reference id and position.
Looking quickly into the source code this is controlled by the -w arguments.
It defaults to 20, which is why you don't get an output for your 24 matches.

It means the read is unmapped. No origin in the reference could be found that was sufficiently similar to the read to call it a proper alignment.
In the SAM specification, * is always used if the information for that field is not available. Therefore, if no mapping information is available (=unmapped) one sets *.


text:
Making some columns if a dataframe numeric
I have two dataframes
> head(a[1:4,1:4])
   Tumor_Sample_Barcode Chromosome Position End_Position
1:                   A1       chr4 90169866     90169866
2:                   A1      chr11 60235747     60235747
3:                   A1       chr1   983023       983023
4:                   A1       chr1 11346060     11346060
>

> head(c)
   Chromosome Position      VAF
1:       chrM     6691 0.610284
2:       chrM    14503 0.693325
3:       chr1 31412236 0.645161
4:       chr1 55693305 0.602941
5:       chr1 69963412 0.709302
6:       chr1 72608266 0.720000
>

> str(c)
Classes ‘data.table’ and 'data.frame':  44175 obs. of  3 variables:
 $ Chromosome: chr  "chrM" "chrM" "chr1" "chr1" ...
 $ Position  : num  6691 14503 31412236 55693305 69963412 ...
 $ VAF       : num  0.61 0.693 0.645 0.603 0.709 ...
 - attr(*, ".internal.selfref")=<externalptr>

I want to merge them like this but complaining about not being integer
> a=a[c, on = c("Chromosome","Position"), VAF := i.VAF]
Error in bmerge(i, x, leftcols, rightcols, roll, rollends, nomatch, mult,  :
  Incompatible join types: x.Position (character) and i.Position (integer)

I have tried something like these but not working
> c=c[!is.na(as.numeric(as.character(c$Position))), ]
>
> c=c[!is.na(as.numeric(as.character(c$VAF))), ]
>
> c$VAF <- as.numeric(as.character(c$VAF))

c$Position <- as.numeric(as.character(c$Position))

which not working

Please add output of str(head(a[1:4,1:4]))

This type error is the same you had here : Merging files gives memory allocation error. As I wrote in comment, updating the type of one data.table column is done this way :  c[, position := as.integer(position)]

Your issue seems to be related to the column Position being a character column in a, you should try:
library(data.table)
a$Position <- as.numeric(a$Position)
a[c, on = c("Chromosome","Position"), VAF := i.VAF]


text:
Seurat FindMarkers() output, percentage
I have generated a list of canonical markers for cluster 0 using the following command:
cluster0_canonical <- FindMarkers(project, ident.1=0,  ident.2=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14), grouping.var = "status",
min.pct = 0.25, print.bar = FALSE)

And I got the following output, but I do not understand why the canonically expressed ribosomal markers seem to be expressed at high levels in both the cells within the cluster AND cells that are NOT in the cluster (i.e. both pct.1 and pct.2 are close to 1, but P value is highly significant). We would not expect this for canonical markers - we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster. I don't seem to have such an issue with non-ribosomal genes within cluster 0 and cluster 0 is the only cluster where I get a high number of ribosomal genes as canonical markers.

Your interpretation of pct.1 and pct.2 is misguided. If pct.1 == 1, the given gene was detected in all of the cells. This carries no information about high expression or low expression in the cell, just that it was detected. avg_logFC is the indication of expression levels. So your expectation, that

"we would instead expect the genes to be expressed at high levels in the cells within the cluster but at low levels in cells not within the cluster"

holds true, as indicated by the positive avg_logFC and significant p-values of the listed RP genes.
Seurat use nautral log, so the FC of RPS6 in cluster 0 vs. all other clusters indicated is 2.718281828459^.55947=1.750.

Seurat FindMarkers() documentation.


text:
Existing tool for converting gff3 to genbank (gbk)
I want to convert my gff3 annotation files to genbank format for use in Mauve. I found the seqret tool here https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/ which can perform this task, but my files (bacterial genomes) are too big.
Can anyone suggest a simple method to make this conversion?

Hi there. I don't know any simple method for what you ask, mainly because gff3 format does not contain sequences (only annotations) while fasta & genbank do. This thread on biostars might help tho :  Gff3 + Fasta To Genbank. As the title suggest, you will need your genome gff3 and fasta reference as inputs.

THanks for taking the time to consider my question Paul.  My gff3 files have the sequence in FASTA at the end of the file (following a ##FASTA indicator), but i take your point - will need some further leg work on this one!  Thanks for the link also.

I believe seqret is already the simplest approach. If you install the emboss command line utility, genome (file) sizes will not be a problem. You can even install it with conda, see https://anaconda.org/bioconda/emboss.
After that, you can use it as the example:
seqret -sequence {genome file} -feature -fformat gff -fopenfile {gff file} -osformat genbank -osname_outseq {output prefix} -ofdirectory_outseq gbk_file -auto

Hope it helps

You can use BioPerl to do this. Read in the GFF and make gene objects and write out those gene objects in Genbank format.
Perhaps this tool (doing the opposite) will get you started:
https://metacpan.org/pod/bp_genbank2gff3.pl
Here is a script doing what you ask but it doesn't use any libraries (so is probably buggy):
https://gist.github.com/avrilcoghlan/5386810
Hope that helps,
Dan.

Here's a python script based on the bcbio-gff (reading gff3) and biopython (writing genbank): https://github.com/chapmanb/bcbb/blob/master/gff/Scripts/gff/gff_to_genbank.py

just random ref, can pip install bcbio-gff and pip install BioPython to get these dependencies for this


text:
Replace a prefix of a string
I have a data frame like this
> head(a[1:2,1:4])
   Tumor_Sample_Barcode Chromosome Position End_Position
1:                   A1       chr4 90169866     90169866
2:                   A1      chr11 60235747     60235747
    > dim(a)
[1] 753655    134
>
> unique(a$Tumor_Sample_Barcode)
 [1] "A1"  "A10" "A11" "A12" "A13" "A14" "A15" "A16" "A17" "A18" "A19" "A2"  "A20" "A21" "A22" "A23" "A24" "A25" "A26" "A27" "A28"
[22] "A29" "A3"  "A30" "A31" "A32" "A33" "A4"  "A5"  "A6"  "A7"  "A8"  "A9"

I want to change a$Tumor_Sample_Barcode  from A1:A21 to res1:res21 but I don't how to do that

I think you are looking for gsub function:
a$Tumor_Sample_Barcode <- gsub("A","res",a$Tumor_Sample_Barcode)

If you prefer, you can use the str_replace_all function from stringr package. But it tooks a little bit more time than gsub:
Dummy data
seq = as.factor(paste0("A",1:21))
library(data.table)
dt <- data.table(Barcode = rep(seq,10000),
                 Value = rnorm(210000))

Test of gsub vs str_replace_all
> system.time(dt$Barcode <- gsub("A","res",dt$Barcode))
   user  system elapsed
  0.077   0.000   0.077

> system.time(dt$Barcode <- str_replace_all(dt$Barcode,"A","res"))
   user  system elapsed
  0.123   0.004   0.126
```

Other options:

Drop 1st character and prefix with "res"
We get 5x speed up if we set fixed = TRUE for gsub.

microbenchmark::microbenchmark(
  # dc37 answer
  gsub("A","res", dt$Barcode),
  stringr::str_replace_all(dt$Barcode, "A", "res"),
  # other options
  paste0("res", substring(dt$Barcode, 2)),
  gsub("A","res",dt$Barcode, fixed = TRUE),

  # bench set up
  check = "identical",
  unit = "relative"
  )

Unit: relative
                                             expr      min       lq     mean   median       uq      max neval
                     gsub("A", "res", dt$Barcode) 5.171682 5.132816 4.283014 5.033968 3.177772 2.545173   100
 stringr::str_replace_all(dt$Barcode, "A", "res") 6.914598 6.870849 5.747396 6.728126 4.809257 3.518024   100
          paste0("res", substring(dt$Barcode, 2)) 1.983757 1.991876 1.779557 1.968488 1.480892 3.415619   100
       gsub("A", "res", dt$Barcode, fixed = TRUE) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000   100


text:
NA12878 equivalent for somatic variant calling
I'm building a somatic variant calling pipeline with snakemake. The pipeline is hosted on Github, and I'd like to use Travis CI for continuous integration testing. Due to privacy concerns, I am looking for a publicly available tumor-normal dataset similar to Genome in a Bottle's NA12878 dataset that I can place online without privacy concerns and use for integration testing. Any sources to find tumor-normal FASTQ data that is publicly available?

Are you after whole exome sequencing (WES) or whole genome sequencing (WGS) for you variant calling pipeline?

Right now the pipeline only supports WES data, but I'm planning to add WGS support so both data types are of interest!

If the site from the suggested duplicate above doesn't work, do try emailing the study's authors. I did so last year and they were incredibly helpful and actually put the site back up for me! If that isn't possible, leave me a comment here (use @terdon) and I can try and get the data to you.


text:
Module preservation and hub genes finding
I ran module preservation analysis between two condition and on the basis of preservation z summary i would to go for downstream analysis.So I m interested in looking for hub-genes of those modules which have lowest z summary. But so far in the tutorial i didn't find a section in the WGCNA tutorial where how to find hubgenes after module preservation analysis was mentioned.
Any suggestion would be really appreciated

Simply select the module(s) you are interested in and look for hub genes using the standard calculations. Having done a module preservation analysis does not change the procedure in the slightest.

but so far what i used to do i run network analysis on individual condition and find hub genes but in module preservation i haven;t found the step or the output from which I can look for hub-genes it would be really helpful ,if you can cite me "using the standard calculations. " in your site

Well, if you ran module preservation, you must have already previoulsy done a network analysis of at least the reference data. Go back to that network analysis and focus on the modules flagged by the preservation analysis. Use kME in those modules to select the top hubs. Section 3d of Tutorial I shows you how to create a table of kMEs.

thank you for the reference