text:
Tool or script to parse annotated VCF files
I've annotated VCF files using snpEff and looking for a tool or script to parse the VCF file and clean up the file to make it interpretable for a biologist.

(edit)you can filter the VCF annotations with snpsift, I've also written a VcfFilterSequenceOntology http://lindenb.github.io/jvarkit/VcfFilterSequenceOntology.html
I've written vcf2table: http://lindenb.github.io/jvarkit/VcfToTable.html
It decodes VEP and SNPeff annotations:
>>chr1/10001/T (n 1)
 Variant
 +--------+--------------------+
 | Key    | Value              |
 +--------+--------------------+
 | CHROM  | chr1               |
 (....)
 VEP
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 | PolyPhen                 | EXON | SIFT           | ALLELE_NUM | Gene            | SYMBOL | Protein_position | Consequence                                   | Amino_acids | Codons  | Feature         | BIOTYPE              |
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 | probably_damaging(0.956) | 8/9  | deleterious(0) | 1          | ENSG00000102967 | DHODH  | 346/395          | missense_variant                              | R/W         | Cgg/Tgg | ENST00000219240 | protein_coding       |
 |                          | 3/4  |                | 1          | ENSG00000102967 | DHODH  |                  | non_coding_exon_variant&nc_transcript_variant |             |         | ENST00000571392 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000572003 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000573843 | retained_intron      |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  |                  | downstream_gene_variant                       |             |         | ENST00000573922 | processed_transcript |
 |                          |      |                | 1          | ENSG00000102967 | DHODH  | -/193            | intron_variant                                |             |         | ENST00000574309 | protein_coding       |
 | probably_damaging(0.946) | 8/9  | deleterious(0) | 1          | ENSG00000102967 | DHODH  | 344/393          | missense_variant                              | R/W         | Cgg/Tgg | ENST00000572887 | protein_coding       |
 +--------------------------+------+----------------+------------+-----------------+--------+------------------+-----------------------------------------------+-------------+---------+-----------------+----------------------+
 Genotypes
 +---------+------+-------+----+----+-----+---------+
 | Sample  | Type | AD    | DP | GQ | GT  | PL      |
 +---------+------+-------+----+----+-----+---------+
 | M10475  | HET  | 10,2  | 15 | 10 | 0/1 | 25,0,10 |
 | M10478  | HET  | 10,4  | 16 | 5  | 0/1 | 40,0,5  |
 | M10500  | HET  | 10,10 | 21 | 7  | 0/1 | 111,0,7 |
 | M128215 | HET  | 15,5  | 24 | 0  | 0/1 | 49,0,0  |
 +---------+------+-------+----+----+-----+---------+


text:
Removing PCR duplicates in RNA-seq Analysis
After reading some of the forum posts in Biostar and SeqAnswers I find it very confusing whether to filter out the duplicate reads from aligned files or not. As far I understand it's very difficult to distinguish between highly expressed genes and duplicate reads and we may lose important information during the filtration process.
So, is it really necessary to remove the duplicates in differential expression analysis using RNA-seq data?

It all depends on how severe the problem is. In general you should only remove them if it is a very high percentage.

Any particular threshold like more than 35%?

Unfortunatly not. The problem is by removing them you might remove real duplicated data thereby under powering your analysis - but by not removing them you might count PCR artefacts thereby increasing false positive. So either option have a potential bias.
I would probably try doing your analysis both with and without removal and then manually inspect some of the results only found by one approach to se which you find more trustworthy. What you should pay special attention to are whether the coverage is more uniform or seem to be affected by loads of reads from single positions.

For normal RNA-seq PCR duplicates are normally kept in, but the duplication rate can be used as a quality control: The higher the duplication rate, the lower the quality. For expression analysis, it is probably best to discard high duplication rate samples, rather than deduplicate them.
In general, the smaller the amount of RNA input into the library prepartion the worst the duplication. Many protocols for very low input quantities (such as single cell) include random barcodes called UMIs (Unique Molecular Identifiers). These allow PCR duplicates to be distinguished from genuinely independent molecules that just happen to have the sample mapping position.

Since the OP does not mention which RNA classes the analysis, let me add that UMIs are also extremely important for sRNA, in particular piRNAs, which quite often map to the exact same locations in the genome and have the same sequence. As afar I know current methods to mark duplicates in mapped reads would result in extreme under-reporting of these RNAs.

Generally you should just leave them as is. One does remove/mark duplicates in DNA seq.
For further read check this Nature paper

Your answer could be improved by hinting on why, according to the paper you cite, removing "PCR duplicates" shouldn't be done in RNASeq. E.g., by citing the one crucial sentence of the abstract: "We find that a large fraction of computationally identified read duplicates are not PCR duplicates and can be explained by sampling and fragmentation bias."

Sure..I could have excuse me for short answer.

Some of the researchers in my institute have researched on this and use UMIs (Unique Molecular Identifiers) to eliminate PCR duplicates. They have been able to increase the reproducibility of the rna seq analysis they had been working on. The UMIs capture some lengths of diverse nucleotides from the rna-sequence. I cannot validate how much UMIs can help in deduplication, as even Nature journal says that UMIs specifically work on the 3’ end. But it has been observed that eliminating duplicates doesn’t improve the accuracy of quantification and it could introduce a bias in the data. So do not eliminate the duplicates, it's just a waste of time.

Hello Shirin, I think it's important to provide a reference for the following statement in your answer: "But it has been observed that eliminating duplicates doesn’t improve the accuracy of quantification and it could introduce a bias in the data."

source: Stack Exchange

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

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

sort: votes


text:
Where can I find summary data for how common certain mutation *types* are?
I'd like to know how common certain mutation types are in public data sets like the 1000 Genomes, ExAC, and ESP6500. Specifically, I'd like to know the distribution of stop-gains, stop-losses, frameshift, and other mutation types. For example, what is the median number of stop-gains observed across individuals in the ExAC population? I'd like to have the full distribution, but median values would work.
To be clear, I am not looking for allele frequencies.
I didn't find this kind of summary data when I looked at the respective websites. Is there a good source for these summary data?

Are you looking for statistics on the most severe variant at each genomic location? It is possible for many different transcript variants to map to the same location, which can influence the results that come out of this.

Does Figure 4A help? https://www.nature.com/nature/journal/v491/n7422/fig_tab/nature11632_F4.html

You might be able to calculate yourself the stats you want if you download the data.

For some types of variant annotations, you can pull all the SNPs from dbSNP as they do in this paper.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3500177/

An important clarifying question: Do you only care about this information in the context of humans?

This dbsnp github repo seems to offer that information.
https://github.com/ncbi/dbsnp/tree/master/Build%20Announcements/151

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


text:
Basic questions about GSEA
I'm currently learning about Gene Set Enrichment Analysis (GSEA) in the hopes of using it in my analysis of differentially expressed genes, and I just had a few questions about the program, specifically about GSEAPreranked, which I need cleared up.

On the ranked list needed for GSEA input, should the list include all genes, or only those that pass a certain threshold of significance (i.e. fold change higher than 2, p value less than 0.05, etc.)? Ideally I'd like to sort the genes by fold change alone as I don't trust my p values as much, so should I only include genes with high fold changes?
I am comparing multiple conditions of disease with different treatments. Am I correct that GSEA only compares two conditions? If this is the case should I run GSEA for each control/treatment comparison? Would this be conventional?

You seem to refer to the GSEA provided by the Broad institute, (there are other GSEA algorithms).
1) You can provide whatever you wish, but if you want to know if those gene sets in which side of the ordered list are they, then provide all the list (of genes) you have.
2) GSEA analyize if the order of a given list distributes in certain way the elements of another (unordered) list. In your case you would need to compute the GSEA for each treatment vs control comparison.

Use the following R package for Gene Set Enrichment analysis of RNA-seq data: seqGSEA
There is another R package (fgsea) recently published called "Fast Gene Set Enrichment Analysis" by Alexey Sergushichev.

According to the GSEA documentation:

we are recommending use of the GSEAPreranked tool for conducting gene
  set enrichment analysis of data derived from RNA-seq experiments

Then more specifically to your question:

Based on your differential expression analysis, rank your features and
  capture your ranking in an RNK-formatted file. The ranking metric can
  be whatever measure of differential expression you choose from the
  output of your selected DE tool. For example, cuffdiff provides the
  (base 2) log of the fold change.

GSEAPreranked is a limited version of GSEA.
If your sample number is enough(more than 5,6 samples in each condition), you can use GSEA version.
Moreover, you can ask this question to the gsea-help google group.

please add links to support your answer


text:
How can I use Arlequin via the command line?
I've got a decent knowledge of programming (incl. bash scripting) but I fail to understand how Arlequin works. Could you please help me with a very simple reproducible example on how to use Arlequin via the command line?
As the .zip file comes with its pre-compiled version, I'll just note that I'm on Mac OS X but of course, a reproducible example with Linux (or at worst with Windows) would also be helpful.

In the manual, I read in the "systems requirement": "Windows
XP/Vista/7". I guess that what you want to use is the arlumstat thing described here: http://cmpg.unibe.ch/software/arlequin35/man/arlsumstat_readme.txt

But in this readme file, it says you need a configuration file generated by the windows version of arlequin. This software looks rather messy. If you know the formulas to compute your statistics, you might in the long run benefit from implementing the calculations by yourself or find a bioinformatician friendly alternative package.

I don't know if Arlequin is messy but the version 3.5 alone has been cited 15274 times! I sure could implement the technic myself but I specifically want to use some commonly used software. I don't think there exist other implementations of FDist since Lositan "disappeared"!

I guess all these people had a windows version to start with. Searching for "Arlequin linux" got me the following page: https://trapa.cz/en/arlequin-r-linux This person confirms that he had to run the windows version first. He did it using wine, a kind of windows emulator for linux. I don't know what the options are for MacOSX. I will definitely not call this program "bioinformatician-friendly".

@bli ok! I am surprised that Arlequin is not more user friendly. A reproducible example with Windows will still be helpful. Thank you

Note that it may well be user friendly, for users that have windows. That's not the case of all bioinformaticians.

A follow-up question can be found here :)

For MAC OSX, the executable found in arlsumstat_macosx (arlsumstatmac_64bit) seem to work appropriately. The directory for arlsumstat_macosx does not contain any examples but you can use example files from arlecore_macosx/Example files_linux using this syntax.
The complication is that arlsumstatmac_64bit will look in the current directory for the file arl_run.ars. There is a file with the same name in arlecore_macosx but it is not compatible. So you have to be in the arlsumstat_macosx.
In short, cd to arlsumstat_macosx and do
./arlsumstatmac_64bit path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.ars path/to/arlecore_macosx/Example\ files_linux/Freqncy/cohen.arp

, where the .arp file contains the data and the .ars file contains the settings (the description of what you want Arlequin to do).

A concise description for using Arlecore on terminal for MAC, LINUX or Windows is provided here:
http://cmpg.unibe.ch/software/arlequin35/man/arlecore_readme.txt

Hi there, thanks for your contribution the link seems useful, but OP specifically asked for an example. It would be nice to provide one.


text:
what percentage of the human genome is MAPQ=0?
When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, what percentage of the non-N fraction of the human genome is MAPQ=0? This is, what part corresponds to regions that can't be uniquely mapped with 2x150bp reads.
And, how many genes are affected by MAPQ=0 regions?
I presume the numbers will dance around depending on fragment sizes, read quality, etc. but I am happy with some starting numbers.

Are you expecting someone to have these numbers handy or do you want to know how to figure this out? :)

UCSC used to have downloadable mappability tracks. These don’t seem to exist any more; the best bet now is to generate them yourself, using GEM.

You can start with the repeat masked regions of the genome.  That would reduce the search space if you want to look at 300bp regions.  What is the insert size of your 2x150 library?  That will play a large role in the mappability of the pair.

Mapping quality is determined by the repetitiveness of the genome, the sequencing error rate, insert size, the capability of the mapper and the nasty heuristics behind the mapper. MAPQ=0 to one mapper is not necessarily MAPQ=0 to another.
That said, I get what you mean. You want to know the uniqueness/repetitiveness. It is still hard if you want to get a useful answer. For 150bp reads, each reference position is covered by 150 reads. What if 50 of them have no other exact hits elsewhere, but the rest 100 have? Is this position a repeat or not? In addition, what if the 50 are unique only because one mismatch? If there is a variant at that mismatch, the 50 would become repeats or mapped elsewhere.
My preference is to say a position is "unique" under k-long reads if over k/2 reads overlapping the position have no other perfect or 1-mismatch/1-gap hits elsewhere in the genome. Under this definition, 79.3% of human genome are unique for 35bp reads. 92.4% for 75bp reads. I don't have the number for 150bp reads. I guess it will be around 95%. Empirically, 94-95% of human genome is callable with 100bp paired-end reads.
As to other measurements, the most common one is the fraction of reads that has an exact hit elsewhere. You can use Fabio's method. It gives a good enough estimate once you simulate over 1 million reads. This fraction is around 86% for 35bp reads and 95% for 75bp reads, as I remember. The problem with this approach is this fraction is not very informative to variant calling due to the issues I talked about. Another way is to use RepeatMasker. It is worse. RepeatMasker masks 50% of human genome, but excludes segmental duplications where short reads can't be confidently mapped.

"For 150bp reads, each reference position is covered by 150 reads." sounds wrong to me. This is what I think should be in the sentence "For 150bp reads, each read mapped to reference has a 150 positions."

@KamilSJaron think about 150X coverage.

Ah, all right, now I get it. From the sentence, I thought that this supposed to be a property of 150bp reads.

As has been said before, mappability to the 'human genome' depends on a number of factors, among these the reference version and type of reads, for which you are interested in GRCh38 and 2x150bp reads. Although I am not aware of numbers accounting for these particular reference and type of Illumina reads, the 1000 genomes project has provided the community with a similar and close estimation that you might be interested in considering regarding your inquiry.
Similar to your question, the 1K project estimated 'the proportion of the human genome that is less accessible to short reads'. In these estimates the human genome is GRCh37 and the types of reads in question are mostly 2x Illumina with a mixture of lengths with the longest being up to 250bp. In these estimates each base in the human genome is considered (and marked) 'less accesible' according to these criteria:
L - depth of coverage is much lower than average
H - depth of coverage is much higher than average
Z - too many reads with zero mapping quality overlap this position
Q - the average mapping quality at the position is too low
Each of these criteria has "standard" and "strict" thresholds for a base to be considered - or not - in each category. You can read more in the link below:
ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/supporting/accessible_genome_masks/README.accessible_genome_mask.20140520
According to the strict thresholds, the human genome has about 16.8% of "Z" and 3.1% of "Q" bases, respectively. Considering the "Z" and "Q" criteria as a proxy for ~ mapq=0, about 19.9% of the human genome can not be uniquely mapped.

I would suggest you to run some simulation. You can use wgsim to simulate reads that are as much similar to what you want an answer for.
I don't think there are faster methods.

I don’t think this is an appropriate approach to investigate mappability because the sample will be too small. Tools like GEM perform exhaustive calculations.

You might want to check out this paper:

Ebbert MTW, et al., Systematic analysis of dark and camouflaged genes reveals disease-relevant genes hiding in plain sight. Genome Biol. 2019 May 20;20(1):97. doi: 10.1186/s13059-019-1707-2. PMID: 31104630; PMCID: PMC6526621.

They have collected a list of "dark" (regions that have few reads aligning to them) and "camouflaged" (regions that have ambiguous alignment) regions of the human genome:

Here, we identify regions with few mappable reads that we call dark by depth, and others that have ambiguous alignment, called camouflaged.
[. . .]
Other dark regions arise, not because the sequencing is inherently problematic, but because of bioinformatic challenges. Specifically, many dark regions arise from duplicated genomic regions, where confidently aligning short reads to a unique location is not possible; we term these regions as “camouflaged”. These camouflaged regions are generally either large contiguous tandem repeats (e.g., centromeres, telomeres, and other short tandem repeats), or a larger specific DNA region that has been duplicated (e.g., a gene duplication) either in tandem or in a more distal genome region. In fact, many genes in the human genome were duplicated over evolutionary time and are still transcriptionally and translationally active (e.g., heat-shock proteins) [3–9], while others have been duplicated, but are considered inactive (i.e., pseudogenes). Regardless of whether the duplication is active, however, any genomic region that has been nearly identically duplicated and is large enough to prevent sequencing reads from aligning unambiguously will be “dark”, because the aligner cannot determine which genomic region the read originated from.

The authors have made their scripts available on github, along with bed files for the relevant regions. The answer to your question, i.e. a specific percentage, is hard to give since it varies enormously with the genome build and sequencing technology, but this paper should have all the details you want. Have a look at Table 1:


text:
How can I compare two bed files?
I am trying to translate (lift over) bed files describing genomic regions from hg37 to hg38. I have tried both UCSC's LiftOver tool and CrossMap but saw that they give me different results. I therefore need a way of assessing how correct the results are.
I decided to test by downloading a bed file of the RefSeq genes for hg37 and one for hg38 from UCSC. Then, I run liftOver to translate the hg37 file to hg39 coordinates and now I want to compare the results of the liftover to the bed file I downloaded for hg38.
How can I usefully compare these files? I am hoping for a tool that can report how similar the two are. Ideally, I would like to see three things in the output:

Matching regions: those that are shown with identical positions in the liftover results and the UCSC bed file.
Non-matching regions: those whose coordinates do not match between the two files.
(if possible) fuzzy matches: those that are not identical but pretty close (off by a few nucleotides).

I thought it would be easy to do with a little awk script but it gets quite complicated because certain RefSeq genes are mapped to multiple positions (e.g. MIR4454). And the 3rd requirement (fuzzy) is not trivial to do in awk. I can live without that though, if a tool can give me the other two.

From a set operations viewpoint, consider your hg37 regions as a "reference" set and hg38 as a "map" set.
If you have these sets in BED format, BEDOPS lets you do set operations on them with various overlap criteria.
Specifically, you can use BEDOPS bedmap to map hg38 regions to hg37 regions. The bedmap tool lets you assign overlap criteria from as little as one base to as much as full, mutual overlap.
Once you do mapping, you can do operations, like counting the number of regions that meet the overlap criteria between reference and map sets, and calculate percentages.
For example, use --exact to get exactly-matching regions between hg37 and hg38:
$ bedmap --count --exact hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt

You can do fractional ("fuzzy") overlaps, where (for example) at least 50% of an hg38 region has to overlap an hg37 region:
$ bedmap --count --fraction-map 0.5 hg37.bed hg38.bed | awk '{s+=$1}END{print s;}' > fuzzy_0p5_counts.txt

Getting a percentage out of a "fuzzy" result can be complicated by overlaps between hg38 regions that overlap a reference hg37 region. A more complicated procedure can be implemented with the --echo-map option and some command-line work to implement logic to deal with that case more stringently.
To get a straight-up count of no overlaps, build the inverse — or mask — of your hg38 regions, using Kent tools fetchChromSizes and BEDOPS bedops --difference:
$ fetchChromSizes hg38 | awk '{print $1"\t0\t"$2}' | sort-bed - > hg38.bounds.bed
$ bedops --difference hg38.bounds.bed hg38.bed > hg38.masked.bed

The regions in hg38.masked.bed make up all of the genomic space in hg38, except for your lift-over regions.
Apply --count --fraction-ref 1.0 to count the number of reference (hg37) elements which completely overlap these masked hg38 regions:
$ bedmap --count --fraction-ref 1.0 hg37.bed hg38.masked.bed | awk '{s+=$1}END{print s;}' > no_counts.txt

The documentation for bedmap describes these overlap criteria and operations in more detail.

I wrote ParsEval to handle these types of comparisons. ParsEval reports a variety of similarity statistics at both the nucleotide level and feature (exon) level. It doesn't explicitly support "fuzzy" matches, but it does report scaled values between 0.0 and 1.0 for similarity, so you could get a similar feel by definig a similarity cutoff for filtering.
Caveat: ParsEval doesn't handle BED files, only GFF3.

Ah, if ParsEval also handled BED it would surely be the holy grail of annotation comparison. ;-)

@KonradRudolph Ah, if everyone would use BED for visualization (like it was designed) and GFF3 for generic feature annotation (like it was designed), the world would be a better place. ;-)

@KonradRudolph I joking...partially. But yes, your point is well taken. :-)

I was mainly just trying to make a cheap pun based on the tool’s name. But it’s no secret that I find the G[TF]F format horrible.

I was trying the code from @Alex Reynolds, however, I got retrieved an error:
bedmap --count --exact SRR800856_*.bed | awk '{s+=$1}END{print s;}' > exact_counts.txt
May use bedmap --help for more help.

Error: Option SRR800856_ngs_te_mapper_nonredundant.bed does not start with '--'

But when I tried with just two bed files:
bedmap --count --exact SRR800856_ngs_te_mapper_nonredundant.bed SRR800856_popoolationte_nonredundant.bed

I don't get any error.
Could you tell me what I should add.
Thank in advance for your time.
P.D. As I am new here I couldn't get to comment on the specific thread.

The bedmap application only takes two arguments: a reference file, and a map file.
In a certain use case, it can take one file as an argument, which acts as both a reference and map dataset (the reference is mapped against itself). However, that's probably not relevant here.
Basically, you can't pass a wildcard argument to bedmap — for bedmap you must specify either 1 or 2 files.
Please ask a new question to separate what you're trying to do from this parent question (and tag it with bedops) and I should be able to find it, take a look and see if I can help.

Thank you! I will post a new question.


text:
At what processing step should library strandedness type be taken into account?
Suppose I have single-end RNA-seq data for which the reads in the fastq file are reversed with respect to the original extracted RNAs.
Suppose I have the following workflow:

Map the reads on the reference genome, sort and index the resulting sam output.
Use the bam file to make a bigwig file showing reads mapping on the plus and minus strands of the chromosomes.
Use the bam file and an annotation file to quantify the gene expression.

Steps 2 and 3 depend on step 1, but are not dependant on one another.
I'm sure that strandedness needs to be taken into account no later than at step 3.
The questions are the following:

Is it common to have mapping tools with an option to "correct" for strandedness?
Assuming the bam file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Regarding this second point, I'm split between two attitudes:

There should be no correction, the bigwig files are there to represent the "raw" mapping results.
The bigwig files are there to give an idea of the origin and abundance of the originally extracted RNAs, so a strandedness correction should be used if necessary.

Which one is the most common attitude?
If the second one is preferred, then wouldn't it be "safer" / "cleaner" to handle strandedness as early as possible? Possibly at the mapping, or even using a pre-processing step?

I usually take into account library strandness during the mapping to the reference step (but I never worked with bigwig files myself).
Assuming RNA-seq data from Illumina, you can use Hisat2 for alignment to the reference genome, and you have the option --rna-strandness to specify the strand-specificity information (default unstranded). This is particularly important when you wish to have XS strand attribute in your bam file for easier compatibility with downstream software (especially important in my experience if you are planning to use some utility from Cufflinks for your step 3).

Good to know that histat2 has such an option, but if I understand correctly, the aligned reads are not "corrected for strandedness": it is up to the next steps to do this, taking into account the XS attribute. Am I right?

@bli What would "corrected for strandedness" even mean? The aligners just use these options to preferentially produce results in a given manner (i.e., in a given orientation relative to genes).

@DevonRyan I mean that, knowing that the reads are produced by a protocol giving the reverse of the RNAs, the aligner reverse them and report the mapping of the reversed reads. This would amount to pre-processing the reads, and then assume that the reads correspond to the RNAs in the rest of the pipeline. When you say "in a given orientation relative to genes", I suppose that implicitly, an annotation file is given to the mapper, and not just a genome index.

@bli Correct, the strandedness settings are only useful if the aligner is aware of the transcriptome (either due to directly aligning against it or having been given an annotation).

If your aligner has a strandedness option then go ahead and use it. The general idea here being that you'll preferentially align correctly to genes.
Whether to bother with strandedness here depends on your goal. Are you interested in anti-sense transcription? Do you need to use the bigWig track to accurately quantify translational pausing or some other strand-specific phenomenon? If yes, then go ahead and make different bigWig files for each strand. If you just need the bigWig file for convenience in IGV or creating heatmaps (e.g., with deepTools) then you'll generally be fine without accounting for strand. For well-annotated species, the rate of anti-sense transcription is generally low and the rate of genes overlapping isn't absurdly high. Consequently, you can usually get away with using an unstranded bigWig.
If you're using featureCounts you'll certainly want to use the appropriate strandedness setting.

It's relatively common for RNAseq aligners/mappers to offer strandedness options. Certainly tophat2 and hisat(2) do, but also things like salmon.
Again, for bigWig files the question becomes whether you really need two files to get the information you need. You're not going to delete the BAM files, so you can always make them if you need them.

Firstly, right at the start if the experiment it's important that your RNA samples are processed with a strand-specific protocol (e.g. Illumina's TruSeq Stranded) in order to produce stranded libraries for sequencing. If the samples haven't been treated as such, then there's nothing you can do to correct it.
Secondly, I'm not sure stranded protocols are compatible with single-end reads. The aligners I'm aware of use the direction of aligned read pairs to assess the strand from which the RNA fragment has come from. How that works with single-end reads, I don't know?
Finally, to answer your specific question, essentially step 1 is where you need to start with strand-specific options and do so at all subsequent steps as appropriate.

All stranded protocols are compatible with single-end reads. The library prep is identical, you just aren't flipping things and sequencing the other end of the fragment.

@DevonRyan and at the feature counting stage? You can simply trust the reads are aligned to the correct strand, right?

You can trust that reads are aligned correctly given some MAPQ threshold, so if they have an orientation compatible with a given gene then yes you might as well believe that.

Strandedness should be taken into account at the mapping stage.
If you are presented with only a BAM file, ask for the raw reads. If that's not possible, you can recover the reads (as fastq sequences) from the BAM file and/or remap them, but that conversion requires a lot more effort (e.g. most mappers expect the fastq files from each end to be sorted by read name and in the same order).
Regarding your specific questions...

Is it common to have mapping tools with an option to "correct" for strandedness?

Yes, mapping tools can be provided with strand direction information for mapping. The way that this is done depends on the program, for example bowtie2 has options --fr, --rf, and --ff to take strandedness into account. I haven't done any of my own comparisons to determine whether these options actually make a difference in mapping for Bowtie2. I know that Salmon (which does mapping and quantification at the same time) cares a lot more about strandedness, and will complain/warn if a stranded library type is discovered and not specified.

Assuming the BAM file has been built ignoring library strandedness information, is it good practice to correct for strandedness when building bigwig files?

Yes, but it'd be better to remap with [the correct] strandedness information included. When I generate bigwig files to demonstrate strandedness (even assuming that strandedness is incorporated into the mapping), I still need to do a fair amount of legwork to convert that into bigwig files. The data needs to first be split into four groups (e.g. with samtools view -F 0xXX -f 0xXX):

First read, mapped in the reverse orientation
Second read, mapped in the reverse orientation
First read, mapped in the forward orientation
Second read, mapped in the forward orientation

I then combine two of those groups, depending on the library strandedness scheme. For a typical stranded Illumina run, groups 1 and 4 are combined (e.g. with samtools merge), as are groups 2 and 3. These combined groups are turned into separate BigWig files, representing the "genome-oriented" and "genome-reversed" read pairs.


text:
combining the use of workdir and option --jobscript in snakemake
It seems that in snakemake the script specified after the --jobscript cannot be used properly if a workdir: is specified in the snakefile.
The path of the script specified becomes relative to the workdir defined in the snakefile instead of being relative to the current working directory. I do not know if this is a feature or a bug !!
Nevertheless, it becomes quite painful to workaround this problem. The only way I found was to parse the command line in the snakefile, for example:

snakemake -c "qsub" -j 30 --js ./sge.sh --latency-wait 30 -rp

in order to copy ./sge.sh in the output directory defined by the keyword workdir in the snakefile.
As a side effect, if you specify option for qsub on the command line, for example:

snakemake -c "qsub -e ./logs/ -o ./logs/" -j 30 --js ./sge.sh --latency-wait 30 -rp

the logs folder must be created under the workdir directory.
Is there something I don't understand with the --jobscript option? Am I not using it the right way? Am I going against the best practice?

To quote from the snakemake documentation

All paths in the snakefile are interpreted relative to the directory snakemake is executed in. This behaviour can be overridden by specifying a workdir in the snakefile:
workdir: "path/to/workdir"
Usually, it is preferred to only set the working directory via the command line, because above directive limits the portability of Snakemake workflows.

It's not explicitly stated, but it's sort of implied that all relative paths then become relative to the working directory. I would expect that specifying an absolute path would get around this.
As an aside, in my mind setting the working directory is usually only needed on clusters without a shared file system (presumably shared between the worker nodes, but not with the head node), since there you can't cd $WORKDIR before running snakemake. This is then normally done with your scheduler, in such cases.

thx for your answer. I had seen this statement in the doc. I understand that all path in the snakefile are interpreted relative to the workdir (very practical by the way) but I feel it is still strange behavior that an option of the command line is interpreted the same way. Moreover, setting a workdir simplifies a lot the writing of the inputs and outputs (no need to prefix all files with the output dir you want)

@EricC. I agree in this case, I would also have preferred that the working dir not apply to paths specified on the command line.

yes maybe Johannes Köster will one day change this behavior, except if for some reason I don't get, it has other drawbacks.

It's just a guess, but I think you might want to use --keep-target-files flag to not to adjust paths of targets. From snakemake manual :

–keep-target-files
Do not adjust the paths of given target files relative to the working
  directory.
Default: False

So if I understand right Devon's answer, combination of this flag and setting of workdir could force compute target in different directory that the snakemake is executed in.


text:
changing blast parameters using NCBIWWW module
I have found a blog post with a script that I would like to use for my current research project: link
The script is incredibly fast and produces a smooth conservation plot. In the blog post, the author mentions that it would be totally possible to change the remote blast function to a local database. Instead of making a local database, I have opted to try to change the function to still use a remote blast, but with different parameters.
I found the following resource for learning about the NCBIWWW module: link
Using this, I rewrote the line in the script:
blast_handle = NCBIWWW.qblast(blast_method, "nr", search_gi)

to:
blast_handle = NCBIWWW.qblast(program="blastp",
                                      database="refseq_protein",
                                      sequence=search_gi,
                                      entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ",
                                      expect=20000.0,
                                      alignments=10000,
                                      descriptions=10000)

Strangely, I get the exact same graph as an output even with the different parameters. Am I doing something about this search incorrectly?
Thanks.

Update:
I got my BLAST to match the webpage. my current parameters are:
blast_handle = NCBIWWW.qblast(program="blastp", database="refseq_protein", sequence="P02686.3", word_size=6, expect=20000.0, hitlist_size=10000, gapcosts="11 1", matrix_name="BLOSUM62", filter="F", genetic_code=1, threshold=21, composition_based_statistics=2, alignments=10000, descriptions=10000, entrez_query=" txid2 [ORGN] OR txid4751 [ORGN] ")

With this new blast_handle, the script no longer works. I get the following error codes:

line 191, in
main(sys.argv[1])
line 27, in main
blast_rec = ncbi_manager.remote_blast(protein_gi, "blastp")
line 184, in remote_blast
return rec_it.next()
line 575, in parse
raise ValueError("Your XML file was empty")
ValueError: Your XML file was empty

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


text:
Does Prokka or gff3 change the sequence / subject ID of fasta files?
I used prokka to create a database for some fasta files but I noticed a strange difference between the prokka and fasta files. Normally the seqid for the gff3 output for prokka has seqid's that match the sequence id's of the fasta files. The files I just produced with prokka have a different seqid than my fasta files. My fasta files look like this
>NODE_108_length_645_cov_0.679537_ID_215

However the prokka output gff3 file for the same fasta file uses a seqid like this
gnl|X|JHIOGLGG_28

What would cause this difference? And do the last digits still correlate to the same seqid's even though the string before it has changed?

Prokka is an annotation tool. It takes your contigs file with fasta headers >NODE_108...and searchs for ORFs and their putative functions. Prokka produces submission ready output, the new fasta headers you see >JHIOGLGG_28...are locus tags and gene description from the annotation and are not related to contig node data. Please see the github information page for more info.


text:
What does the "X_OS_IND" column mean?
I want to survival analysis using the subset of TCGA LUAD dataset, which identifier are located here.
> head(LUAD_clinicalMatrix[,c("X_EVENT","X_OS_IND")])
  X_EVENT X_OS_IND
1      NA       NA
2       0        0
3       1        1
4       0        0
5       0        0
6       0        0

> all(LUAD_clinicalMatrix$X_EVENT==LUAD_clinicalMatrix$X_OS_UNIT)
[1] FALSE

I want to know which identifier means the sample is alive or dead, and "1" means alive or dead?

Cross-posted on biostars

According to this link, 1=death 0=censor and null=no data.


text:
Where can I find a GWAS data set about genotypes and phenotypes for Zebrafish and rice?
My goal is to work on GWAS method development as a research project. However, human genome data are usually confidential because of the identification problem, so it's very hard to get them. (application takes a long time and I may not get it at the end.)
Therefore, maybe I can try to play with the data of other species first. So I wonder where I can find some GWAS data for Zebrafish and rice. (The reason I am asking for these two is that I find researchers tend to perform genome analysis to mice, Arabidopsis Thaliana, zebrafish, and rice, and I already found where to download the GWAS data for the first two.)
I can find several data sets (like ZFIN) that have the annotated genome information, but that's only one genome. Where can I download the raw data (like genotypes of 500 zebrafishes and corresponding phenotypes) to start my own analysis?

Possible duplicate of Where can I find a database that has phenotype information together with associated SNPs?

@KamilSJaron Thanks. In my defense, the previous question is asking for the associated results that people already find with their own research. This question, I'm looking for a data set that I can research on. What are the differences between human and these species? I did this mainly because it take a long time to apply for human genome data and we may not get it.

You do not need to defend yourself, I am trying to help you in making question clearer. I would recommend you to [edit] your question to add there more details (including those you wrote in comment). Tell us what you want, what you have tried and why you think it did not work.
If the problem with human is dataset size, why not just subsample human dataset? However, GWAS is notoriously known for requirement of huge sample sizes.

I think it will be more productive for you if you choose one of the public human datasets and practise your GWAS methods on that. A couple of examples examples:

1000 genomes -- use "population" as your trait of association. You'll find that you get association across the whole genome, so the analysis method would change from "what bits are associated" to "what are the bits that are most consistently associated"
NCI Genomic Data Commons -- open access (and controlled access) data associated with cancer

I think you need bit more precise phenotype information than population, if you want to test a potential of method. For population you can just look as 2. and 3. axis of PCR and it is going to work (you do not need any formal statistics)...

But that's not going to tell you what markers are associated. I've had success in the past doing a GWAS-like analysis on populations.

You can look at the projection of loci to principal components and therefore you see what loci are actually responsible for the separation.
I do not doubt that you can associate loci with populations. I am trying to say that you do not need any sophisticated method for that because the signal is very strong. I had an impression that the challenge of GWAS studies lays in identification of weak signals, therefore for benchmarking taking "easy" data is not really useful (maybe for sanity tests).


text:
Counts obtained by featureCounts seem much less than observed coverage
I have surprisingly low counts when running featureCounts on some (single-end) RNA-seq data mapped on C. elegans genome using hisat2.
To more easily show the problem, I generated a small subset of the bam file and of the annotation file I'm using. Here is what I can see when loading these two files on IGV:

The coverage for the "his-11" and "his-15" genes peak at around 3600. For "his-12" and "his-16", it peaks at more than 1000
Here are the content of my annotation file:
$ cat his_11-16.gtf
II  ensembl transcript  13817759    13818143    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001890"; gene_name "his-16"; gene_source "ensembl"; gene_version "1"; p_id "P8354"; transcript_biotype "protein_coding"; transcript_id "ZK131.10"; transcript_name "ZK131.10"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS11539";
II  ensembl transcript  13818371    13818739    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001889"; gene_name "his-15"; gene_source "ensembl"; gene_version "1"; p_id "P8185"; transcript_biotype "protein_coding"; transcript_id "ZK131.9"; transcript_name "ZK131.9"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52436";
II  ensembl transcript  13819626    13819937    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001888"; gene_name "his-14"; gene_source "ensembl"; gene_version "1"; p_id "P21890"; transcript_biotype "protein_coding"; transcript_id "ZK131.8"; transcript_name "ZK131.8"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS50308";
II  ensembl transcript  13820210    13820620    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001887"; gene_name "his-13"; gene_source "ensembl"; gene_version "1"; p_id "P2149"; transcript_biotype "protein_coding"; transcript_id "ZK131.7"; transcript_name "ZK131.7"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS35106";
II  ensembl transcript  13821198    13821581    .   -   .   gene_biotype "protein_coding"; gene_id "WBGene00001886"; gene_name "his-12"; gene_source "ensembl"; gene_version "1"; p_id "P26082"; transcript_biotype "protein_coding"; transcript_id "ZK131.6"; transcript_name "ZK131.6"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS23693";
II  ensembl transcript  13821778    13822314    .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00001885"; gene_name "his-11"; gene_source "ensembl"; gene_version "1"; p_id "P10535"; transcript_biotype "protein_coding"; transcript_id "ZK131.5"; transcript_name "ZK131.5"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS52792";

The featureCounts run:
$ featureCounts -a his_11-16.gtf -o his_11-16_counts.txt -t transcript -g gene_name -O his_11-16_sorted.bam

        ==========     _____ _    _ ____  _____  ______          _____
        =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \
          =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
            ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
              ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
        ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
      v1.5.2

//========================== featureCounts setting ===========================\\
||                                                                            ||
||             Input files : 1 BAM file                                       ||
||                           S his_11-16_sorted.bam                           ||
||                                                                            ||
||             Output file : his_11-16_counts.txt                             ||
||                 Summary : his_11-16_counts.txt.summary                     ||
||              Annotation : his_11-16.gtf (GTF)                              ||
||      Dir for temp files : ./                                               ||
||                                                                            ||
||                 Threads : 1                                                ||
||                   Level : meta-feature level                               ||
||              Paired-end : no                                               ||
||         Strand specific : no                                               ||
||      Multimapping reads : not counted                                      ||
|| Multi-overlapping reads : counted                                          ||
||   Min overlapping bases : 1                                                ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================================= Running ==================================\\
||                                                                            ||
|| Load annotation file his_11-16.gtf ...                                     ||
||    Features : 6                                                            ||
||    Meta-features : 6                                                       ||
||    Chromosomes/contigs : 1                                                 ||
||                                                                            ||
|| Process BAM file his_11-16_sorted.bam...                                   ||
||    Single-end reads are included.                                          ||
||    Assign reads to features...                                             ||
||    Total reads : 35037                                                     ||
||    Successfully assigned reads : 1849 (5.3%)                               ||
||    Running time : 0.00 minutes                                             ||
||                                                                            ||
||                         Read assignment finished.                          ||
||                                                                            ||
|| Summary of counting results can be found in file "his_11-16_counts.txt"    ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

And the resulting counts file:
$ cat his_11-16_counts.txt
# Program:featureCounts v1.5.2; Command:"featureCounts" "-a" "his_11-16.gtf" "-o" "his_11-16_counts.txt" "-t" "transcript" "-g" "gene_name" "-O" "his_11-16_sorted.bam"
Geneid  Chr Start   End Strand  Length  his_11-16_sorted.bam
his-16  II  13817759    13818143    -   385 869
his-15  II  13818371    13818739    +   369 3
his-14  II  13819626    13819937    -   312 953
his-13  II  13820210    13820620    +   411 23
his-12  II  13821198    13821581    -   384 423
his-11  II  13821778    13822314    +   537 1

Why such a low assignment rate, why so little counts?
I see that hisat2 decided to split an important proportion of the reads between orthologs, like "his-16" and "his-12", but my understanding is that the default values for the minOverlap and fracOverlap parameters should ensure that such split reads will be counted:
  --minOverlap <int>  Minimum number of overlapping bases in a read that is
                      required for read assignment. 1 by default. Number of
                      overlapping bases is counted from both reads if paired
                      end. If a negative value is provided, then a gap of up
                      to specified size will be allowed between read and the
                      feature that the read is assigned to.

  --fracOverlap <float> Minimum fraction of overlapping bases in a read that is
                      required for read assignment. Value should be within range
                      [0,1]. 0 by default. Number of overlapping bases is
                      counted from both reads if paired end. Both this option
                      and '--minOverlap' option need to be satisfied for read
                      assignment.

What did I get wrong?

How many of those are labeled as multimappers? Those are excluded by default.

About 62% of the reads appear to have the "secondary alignment" flag bit, and almost 95% have the value for the NH tag higher than 1.

well there's the cause.

From the screenshot it seems that his-11 and his-15 have virtually identical (not merely similar) peak shapes. That should give pause.

@KonradRudolph To me, this indicates that the mapper does a "good" random assignment of reads when multiple equally-scoring mapping locations are possible (his-15 and his-11 have identical sequences). I'm more worried about the reads that are split between two his-* genes. This is most likely wrong.

Given the high level of multimapping in this region, you'll need to use the -M --primary options if you want to keep many of the alignments. I would be very hesitant to use these numbers as input for DESeq2 or similar programs, since it's fairly questionable whether one should fully trust the "randomness" of the aligner's assignments. I'm more comfortable using something like salmon for such cases.


text:
ribo interval list for picard ensembl mm9 gtf
Could some one please help me with understanding how to generate the ribosomal interval list that is required to when using picard metrics this step
java -jar picard.jar CollectRnaSeqMetrics \
      I=input.bam \
      O=output.RNA_Metrics \
      REF_FLAT=ref_flat.txt \
      STRAND=SECOND_READ_TRANSCRIPTION_STRAND \
      RIBOSOMAL_INTERVALS=ribosomal.interval_list

I have a couple of bam files that were aligned using mm9 ensembl gtf as annotation, i did a lot of search and found this link Ribosomal Intervals For Collectrnaseqmetrics but this link or any other source does not explain how to produce or find ribo interval list from ensembl annotation

At least with Ensembl annotations, you can grep rRNA foo.gtf > rRNA.gtf to get the annotated rRNAs. Note however, that the 45S cassette isn't in the mouse reference genome, so you'll never have alignments to 18S, 5.8S, or 28S rRNA. I always get the feeling that the Collect*Metrics programs are made with human sequences in mind, since the cassette is included in the human reference genome.
Anyway, while you can use the annotation and reference sequence to make picard happy, make sure to append the Rn45S sequence to your genome if you want to actually measure rRNA levels. In my experience this is only useful for rRNA-depleted samples, since the rates can vary wildly there.
Edit: To make a file compatible with picard's "interval list", one can do the following with the GTF file:
grep rRNA foo.gtf | awk 'BEGIN{OFS="\t"}{print $1, $4-1, $5}' > foo.bed

According to the picard documentation a BED file can be used as input and the above awk code extracts the appropriate columns and switches to 0-based coordinates used in the BED format.

Interesting, can I directly use the output of grep rRNA foo.gtf > rRNA.gtf as input file? or is there any further modification necessary? I did not understand when you say make picard happy I am using mm9 reference and yes the sequencing done was using ribo depleted technique (could you please point towards some literature or steps)

@novicebioinforesearcher I've updated my reply with the conversion to BED, which allegedly should work. "make picard happy" is another way of saying "give picard a file it can use as input". The results would still not be terribly useful without adding in the Rn45S sequence.

I never came across Rn45S concept (thanks for pointing it out) and have not seen any literature in the tools that use such list what would be your thoughts as to why they would not be speak about it?

@novicebioinforesearcher Many reference genomes (in particular, human and fly) contain it already, so there'd be no need to use Rn45S there. I don't know why it was excluded from the mouse genome, it's one of the things that's always bothered me.


text:
Working with the additional cell lines in NCI-60
NCI-60 has additional cell lines, but when I get results from CellMiner they are not included in there.
Namely, the genes such as:  DLD-1(Colon), KM20L2    (Colon), M19-MEL (Melanoma), don't show in the results of CellMiner.
Is there away to include these additional lines in the correlation calculations?


text:
Influence of pH on protein peptide mass
I'm using the Peptides package in R to calculate the mass of peptides with the mw() method. The thing is that I want to calculate them for different pH values, and thus I wrote this little function to take the pH value into account according this image:

My code is this:
phmw = function(pept_seq , unimod_location, pH ){

  monoisotopic_hydrogen_mass = 1.00794

  if(pH < 7){

    # mol_weight = mol_weight_of_peptide + mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) + monoisotopic_hydrogen_mass

  }else if(pH > 7){

    # mol_weight = mol_weight_of_peptide - mol_weight_of_Hydrogen
    mol_weight = mw(pept_seq, monoisotopic = T) - monoisotopic_hydrogen_mass

  }else if(pH == 7){

    # mol_weight = mol_weight_of_peptide
    mol_weight = mw(pept_seq, monoisotopic = T)
  }

  return(mol_weight)

}

Do you think that this approach is correct or not?

The code is correct as far as I can tell but there’s a lot to be said about style. The code can be made a lot shorter and more readable. Furthermore, don’t use T instead of TRUE, and don’t use return` (which is unnecessary here). Generally I suggest reading the Tidyverse style guide, even if you decide not to follow it.

I'm curious as to the practical application? Is your interest specifically in peptides in solution at different pH values? Just wondering because most peptide MW prediction is with reference to mass spectrometry, where pH is not really relevant.

If you want completely accurate monoisotopic mass you should use the mass of the proton rather than hydrogen. (The difference is the mass of electron).
I second the question of practical applicability since each amino acid will have individual pKa, moreover, in the peptide neighboring amino acids will interact. Thus the protonation and deprotonation will happen at the different pH.

Actually I'm trying to mimic the mass spec for some peptides and see if I can characterize them as detectable by their m/z.
@neilfws Didn't understand exactly what you said. The algorithm behind the mw() assigns the masses of the a.a's from mass spec knowledge ? But one can run a mass spec experiment with different pH values ( Something that changes the mass of the total peptide). Which of these masses does the mw() algorithm uses ?

I'm just wondering whether protons gained by protonation due to low pH are relevant to the ions generated during MS. For example in electrospray MS "it became evident that there is no correlation between solution charging and electrospray charging" - https://www.hindawi.com/journals/ijac/2012/282574/. I would think solution protons would be lost during fragmentation and ion generation.

Your approach in general seems to be OK, but I have a few suggestions:

You are basing all your calculations on whether the pH is above or below 7. However the N and C termini have pKa values that determine whether these groups will be ionized (8 and 3, respectively). You should check whether the pH is above/below these values instead of 7.
Different amino acid side chains can also be ionized and should be accounted for when calculating the MW. Amino acid pKas can be found here
Code wise: you pass a variable unimod_location in your function but it's never used. If it's not being used then take it out.

Amino-acid content (eg Lys, Asp) can have drastic effect on observance in a particular MS-mode.
You should consider including Peptides::pI() or Peptides::charge() to determine the likelihood of (de-)protonation in your calculations.


text:
Solutions for managing data in a small bioinformatics / 'omics lab?
A different sort of problem: even a small 'omics lab generates a lot of data, raw, intermediate and processed. What (software) solutions exist for managing this data, such that "old" data can be retrieved and checked or re-analysed, even after people have left the lab? Important points would be:

ease of installation
ease of putting data in, in an appropriately tagged / labeled fashion (it's no good if the repository is just a centralised bad mess)
useful search and exploration
security (i.e. restricted to members of the lab)
previews / summaries of the data
can accommodate any dataset
local, not SaaS

It seems from my explorations of the topic that there's very little between the primitive (e.g. handrolled Access databases) and big industrial solutions. Things I have looked at include:

Dataverse: very popular, installation seems complex and unclear if uploading is that easy
DSpace: mostly for publications and documents
CKAN
OSF: used this for a while, integrates with a lot of services but uploading data seemed awkward

I don't know of any prebuilt products, but I can describe how we managed this in my postdoc lab, and how I plan to manage it in my newly-started group.
Rule 1: All work happens in the projects project directory on the central
filestore, not on your desktop or laptop.
Rule 2: Heavy computational work is done by the groups standard analysis pipelines. Interpretation is done in jupyter notebooks or Rmarkdown.
A project has a directory on the group's filestore. That directory has a fixed structure:
proj001----raw_data             *
       |
        --external_datasets
       |
        --src                  *
       |
        --notebooks            *
       |
        --pipeline1
       |
        --pipeline2
       |
        -- etc
 * link to a seperate, backed up filesystem.

Pipelines are where heavy duty analysis happens and everyone uses the same standard pipelines, producing the same standard output files and database structure.
So common pipelines might be the mapping pipeline, the readqc pipeline, the differential expression pipeline, the exome pipeline, the chip-seq pipeline etc...
Pipelines have three important outputs: an automatically generated html report, an sqlite database and files in the export directory.
When we archive a project, this is what we save, along with the pipeline's configuration file and log file.
So if I know that Jane Bloggs did an RNAseq 5 years ago on a cell type i'm interested in, if I konw that that was project 5, I know that in the project 5 directory there will be a diff_expression_pipeline directory, and that it will contain an sqlite database called csvdb and that that will have a table called refcoding_deseq_gene_diff and that table will follow a known format. There will be bigwigs in the export directory. Or the BAMs will be in the mapping_pipeline directory.
Of course the problem remains knowing that Jane Bloggs did this RNAseq and what the project was called. We use a Wiki for this, but its not ideal.

We have been organizing data on a "per project" basis, using a GitHub repository for each project. Each project ends up being a paper (written in R Markdown), so you could have multiple "projects" based on the same datasets. See this as an example:
https://github.com/SchlossLab/new_project
The analyses are all managed with a Makefile that goes from loading the sequence data from the SRA, to performing analyses, to rendering the final manuscript. This way it is always clear how things were done. Here is an example of what that can look like:
https://github.com/SchlossLab/new_project/blob/master/Makefile
As far as organizing the raw data and whatnot, I have just used a simple MySQL database hosted on MAMP (Mac version, Linux is LAMP, Windows is WAMP). It has a GUI and everything, and is pretty easy to use. Here is a resource for getting started on that:
http://microbiology.github.io/PDFs/MySQL-Intro-Workshop-Hannigan.pdf
All of this stuff is free to use, which is awesome! Hope this helps! :)


text:
How to convert a binary matrix of gene presence or absence into a fasta sequence
I have a binary matrix of gene presence or absence which looks like: [roary output]
Gene    sample1 sample2 sample3 sample4
fliI        1       1       1        1
patB_1      1       1       1        1
pgpA        1       1       1        1
osmB        1       1       1        1
cspA        1       0       1        1

How can I convert to fasta so it looks like and technically aligned:
>sample1
11111
>sample2
11110

Is it possible to do this in R?

What do you mean by "technically aligned"?

Well technically aligned as in similar to DNA alignments.

It's not a fasta file, but:
> m
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> # Collapse to labeled strings
> blah = apply(m, 2, function(x) paste(x, collapse=''))
> blah
sample1 sample2 sample3 sample4
"11111" "11110" "11111" "11111"
> # Write that to a file with the appropriate format
> cat(paste(mapply(function(x, y) sprintf(">%s\n%s\n", x, y), names(blah), blah), collapse=""), file="some file.fa")

Make sure to change some file.fa. I have no idea why you would want to do all of this, but this will give you the output you want:
$ cat blah.txt
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

Thanks Devon. This looks great. I want it for phylogenetic analysis so its important the layout comes out as in the matrix so its technically aligned.

A mixture of apply, paste and cat will do this:
> data.mat
       sample1 sample2 sample3 sample4
fliI         1       1       1       1
patB_1       1       1       1       1
pgpA         1       1       1       1
osmB         1       1       1       1
cspA         1       0       1       1
> data.conv <- apply(data.mat,2,paste,collapse="") # aggregate along 2nd dim
> cat(paste0(">", names(data.conv), "\n", data.conv), sep="\n")
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

If writing to a file is required, just add a file parameter to the cat command:
cat(paste0(">", names(data.conv), "\n", data.conv),
    sep="\n", file="out_fastaLike.txt")

One could also use cut and a pair of Python scripts:
transpose.py
#!/usr/bin/env python

import sys

for c in zip(*(l.split() for l in sys.stdin.readlines() if l.strip())):
    sys.stdout.write("%s\n" % ('\t'.join(c)))

mock_fasta.py
#!/usr/bin/env python

import sys

for e in (l.split() for l in sys.stdin.readlines() if l.strip()):
    sys.stdout.write(">%s\n%s\n" % (e[0], ''.join(e[1:])))

Then:
$ cut -f2- in.mtx | transpose.py | mock_fasta.py > in.fa

Which gives:
$ more in.fa
>sample1
11111
>sample2
11110
>sample3
11111
>sample4
11111

Why anyone would use R for text parsing is beyond me, but kids these days etc.


text:
Visualization tools for 3C/Hi-C long-range interaction data?
I have identified a genomic region (CpG island) that is enriched for a characteristic of interest in the vast majority of my samples.
After viewing it in the UCSC genome browser, I have confirmed that it does not overlap any annotated gene promoters, gene bodies, or retrotransposons, but am convinced that it is functionally interacting with something long range.
Is there a publicly available visualization tool where I can input my genomic region and get a list of topologically associated domains from 3C or other chromatin-chromatin interaction data?

What assembly is your genomic region of interest in?

The assembly I used was hg19/GRCh37

I would also recommend two very recent Hi-C visualization frameworks (with some public data available in both): HiGlass and JuiceBox.

HiGlass has the result of 7 different TAD callers available. Example using hg38 can be found here: http://higlass.io/app/?config=IPCHmdOQR4CDY2sqj5VJHQ
Another example using hg19 is here: http://higlass.io/app/?config=OBa3x6FiQtSC5g8Mj6u0pQ
The TAD calls are from Forcato et al. 2017.

If you have more than a dozen or even hundreds of regions that you want to compare you could check a companion tool of HiGlass called HiPiler. HiPiler lets you arrange, cluster, and pile up to thousands of genomic regions. Take a look at this 5min video introduction to see if the tool could help you to answer your questions.
Disclaimer: I am the creator of HiPiler.

Depending on the genome you want, you may be able to use the chorogeneome navigator, which is from some of my colleagues. This uses a number of public HiC datasets and allows you to view interactions.


text:
insert size pre and post trimming
I have a problem here with my rna seq data:
Sequencing details:
rRNA was removed, followed by cDNA preparation and generation of stranded libraries using the TruSeq Stranded Total RNA Sample Prep Kit. Sequencing performed on the HiSeq2500 platform (Illumina) to generate 2 × 125 bp paired-end reads
Alignment and preprocessing
reads were aligned using tophat2(std aligner in pipeline at core (--library-type fr-firststrand), unfortunately the original unaligned files were purged and i only have access to this aligned file which i converted to fastq using following steps

samtools sort -n sample.bam -o sample_sorted.bam
bedtools bamtofastq -i sample_sorted.bam -fq sample_1.fq -fq2 sample_2.fq (get a lot of mate skipping errors here)
check for these reads for adapters and trim them and again use hisat2 (with --rna-strandness RF)

trimmed reads
Reads were trimmed by passing parameters --adapters adapters.file --adapter-trim-end RIGHT --length-dist --threads 12 --adapter-min-overlap 7 --max-uncalled 250 --min-read-length 25 -- to FLEXBAR version 2.4
adapters for all samples using multiqc

fastqc sequence length pre trimming

fastqc sequence length post trimming

rmats error

Incorrect readLength. sample.bam has a read length of 114, while readLength param is 125

picard error on aligned bam from hisat2)

ProcessExecutor   2 "Not creating insert size PDF as there are duplicated header names: All_Reads"

snapshot of picard insertsize metrics file

Questions

Is my approach flawed at some step?
why do i get different read lengths error for each sample when i process it through rmats when i specify len as 125?
why is no histogram being produced when using picard insertsize metrics?
is normal to see the graph shift to the right after trimming?

insertsize for all samples (data from picard collectinsertsize compiled using multiqc)

It's not an unreasonable approach. rMATs is rather picky about its input, but it seems you've noticed that.
rMATs can't handle trimmed reads. It also can't handle soft-clipped reads. You might as well not trim. You'll get a lower alignment rate, but your only other choice would be to either exclude the trimmed reads or trim everything down to the same length.
You have a mixture of RF and FR alignments and the plotting script can't handle that.
Assuming you mean the read length distribution, it's actually shifted to the left, not the right. You had a single length before and a long left tail after trimming.

How does one get a mixture of RF and FR, when the library was FR ? doesn't that mean the alignment is doubt full as well?

@novicebioinforesearcher That's a very good question and not one that I have a good answer for.


text:
Getting protein FASTA sequence based on keyword with python
I would like to gather proteins FASTA sequence from Entrez with python 2.7. I am looking for any proteins that have the keywords: "terminase" and "large" in their name. So far I got this code:
from Bio import Entrez
Entrez.email = "example@example.org"

searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)
searchResult = Entrez.read(searchResultHandle)
ids = searchResult["IdList"]

handle = Entrez.efetch(db="protein", id=ids, rettype="fasta", retmode="text")
record = handle.read()

out_handle = open('myfasta.fasta', 'w')
out_handle.write(record.rstrip('\n'))

However, it can get me several terminases from various organisms, while I need only terminase form bacteriophages (specificly Viruses [taxid 10239], host: bacteria). I've managed to get the nuccore accession ids from NCBI of the viruses I am intersted in, but I don't know how to combine those two informations. The id file looks like this:
NC_001341
NC_001447
NC_028834
NC_023556
...

Do I need to access every gb file of every ID and search for my desired protein in it?

Welcome to Bioinformatics! Did you tried to add the " [taxid 10239], host: bacteria" in the search in Entrez? From the retrieved fastas you could also check if the id is of interest for you and then omit those who aren't

I tried, but I don't know if my syntax is correct. Should it be: term="terminase large bacteria[organism]"

You the term search would be the same query you set on the NCBI page, that's '("Viruses"[Organism]) AND (host[All Fields] AND ("Bacteria"[Organism] OR "Bacteria Latreille et al. 1825"[Organism] OR bacteria[All Fields]))'

Found what I was looking for. In:
searchResultHandle = Entrez.esearch(db="protein", term="terminase large", retmax=1000)

I've added:
searchterm = "(terminase large subunit AND viruses[Organism]) AND Caudovirales AND refseq[Filter]"
searchResultHandle = Entrez.esearch(db="protein", term=searchterm, retmax=6000)

which norrowed down my searches to the desired viruses. Granted it's not filtered by host, but by a taxonomy group, but it is enough for my work.
Thank you @Llopis for additional help

Please remember to accept this answer. The system won't let you accept your own answer just after posting it, but you should be able to now.


text:
blasting a refseq protein does not show the protein in the result set
Can anyone explain me, why I don't find a specific protein with a blast that was took before from the NCBI refseq database?
Specifically, I was trying to blast the protein with the accession number "NP_420767" and its sequence, respectively, however that protein does not show up in the results. It not only happens when the standard options are chosen, but also, when "Reference proteins (refseq_protein)" as database in the blast options is selected.
I am really puzzled by that.. Shouldn't blast show the initial refseq entry in the results list as it is part of the refseq database and has the same sequence?

Relevant: Duplicate long hits from PSI-BLAST

NP_420767 is represented by the non-redundant refSeq protein WP_010919826, which has the same amino acid sequence. This is not very clearly annotated, but if you scroll down to the sequence in the GenPept entry for NP_420767, you'll see the following:
CONTIG      join(WP_010919826.1:1..799)

See here for more info.

Ah, I understand.. And in the list after clicking on "identical proteins" on the "WP_010919826.1" site, "NP_420767" gets listed.


text:
Filtering bases based on phred qualities with pysam
Is there a way to filter bases in BAM files based on phred quallities through python's pysam ?
I have a code here that

Takes the nucleobases per position from a BAM file using pysam's pileup function
Saves it in ReverseList and ForwardList based on both strands (i.e forward and reverse), I want to reject those bases that have phred quality below 25 to be not stored in the ForwardList and ReverseList lists so that they are not used for further analysis.
 samfile = pysam.Samfile( filename, "rb" )
 ReverseList = [''] * lenref
 ForwardList = [''] * lenref
 for pileupcolumn in samfile.pileup() :
     for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
               continue
         if not pileupread.is_del and not pileupread.is_refskip:
               if pileupread.alignment.is_reverse: #negative
                      ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
               else:
                      ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

 samfile.close()

Where lenref = 16569(length of mitochondrial genome) and filename is the name of BAM file. I want to filter based of phred qualities of bases.

Am I the only one who is missing the point of filtering alignments by phred scores? Phred scores are assigned during base calling, if you want to filter them out, why after mapping but not before? Also variant callers take into account phred scores, having low phred scores should not worsen quality of calling variants.

I am making a tool for analyzing BAM file so I want to filter it again to improve my results.

It is also reported in previous studies.

Welcome to Bioinformatics! This question seem to need more information to get a good answer (Although the existing one answers how to filter as required). Why do you want to filter the bases? What are you trying to achieve?

The PileupRead object as a query_position attribute, which you can use for this:
for pileupcolumn in samfile.pileup() :
    for pileupread in pileupcolumn.pileups:
         if (pileupread.alignment.mapping_quality <= 15):
             continue
         if not pileupread.is_del and not pileupread.is_refskip:
             if pileupread.alignment.query_qualities[pileupread.query_position] < 10:
                 # Skip entries with base phred scores < 10
                 continue
             if pileupread.alignment.is_reverse: #negative
                   ReverseList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]
             else:
                   ForwardList[pileupcolumn.pos] += pileupread.alignment.query_sequence[pileupread.query_position]

Note the 6th line, which implements a filter with a threshold of 10.


text:
Same transcript coordinates in gtf file, different transcript ID
I have a gtf file from Ensembl, and I noticed that several "transcript" annotations have the exact same coordinates. See for instance the third and fourth transcripts ("Y74C9A.2b.1" and "Y74C9A.2b.4") for this gene:
$ grep "WBGene00022276" genes.gtf | awk '$3 == "transcript" {print}'
I   ensembl transcript  10413   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.2"; transcript_name "Y74C9A.2a.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS7921";
I   ensembl transcript  11495   16793   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.2"; transcript_name "Y74C9A.2b.2"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9220";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.1"; transcript_name "Y74C9A.2b.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11499   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.4"; transcript_name "Y74C9A.2b.4"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS9215";
I   ensembl transcript  11505   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.3"; transcript_name "Y74C9A.2a.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS16820";
I   ensembl transcript  11618   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P24752"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2a.1"; transcript_name "Y74C9A.2a.1"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS43845";
I   ensembl transcript  11623   16842   .   +   .   gene_biotype "protein_coding"; gene_id "WBGene00022276"; gene_name "nlp-40"; gene_source "ensembl"; gene_version "1"; p_id "P22572"; transcript_biotype "protein_coding"; transcript_id "Y74C9A.2b.3"; transcript_name "Y74C9A.2b.3"; transcript_source "ensembl"; transcript_version "1"; tss_id "TSS49800";

What does it mean? Do these two "transcripts" actually correspond to an exon, belonging to two different isoforms?

transcript objects cover the co-ordinates from the start of the first exon to the end of the last exon of a transcript (i.e. an isoform). If two different isoforms share the same first and last exons, but have a different set of internal exons, then their transcript entries will be the same, but the set of exon entires associated with each transcript will be different.
For example, consider the two transcripts below:
                      1         2
             12345678901234567890123
transcript_1 |>>>|----|>>|----|>>>>|

transcript_2 |>>|--|>>|------|>>>>>|

These two transcripts would give the following GTF entires
chr1   ensembl   transcript    1   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    1   5   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    10   13   .   +   .   gene_id "gene_1"; gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   exon    18   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_1";
chr1   ensembl   transcript    1   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    1   4   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    7   10   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";
chr1   ensembl   exon    17   23   .   +   .   gene_id "gene_1"; transcript_id "transcript_2";

As Ian explained, these are different transcripts which happen to have the same start and end positions. You have no information on their exonic structure in that file. However, if you look them up at EnsEMBL, you will see:
Transcript Y74C9A.2b.1 :

1   Y74C9A.2b.1.e1  11,499  11,561  -   -   63
    Intron 1-2      11,562  11,617          56
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

This has 4 exons, two of which are non-coding (UTR). The second starts at
Transcript Y74C9A.2b.4

1   Y74C9A.2b.4.e1  11,499  11,557  -   -   59
    Intron 1-2      11,558  11,617          60
2   Y74C9A.2b.1.e2  11,618  11,689  -   -   72
    Intron 2-3      11,690  14,950          3,261
3   Y74C9A.2b.1.e3  14,951  15,160  -   1   210
    Intron 3-4      15,161  16,472          1,312
4   Y74C9A.2a.1.e3  16,473  16,842  1   -   370

Granted, the two transcripts are very simlar, they both have 4 exons, the first 2 of which are non coding (UTR). However, the 1st exon of transcript
Y74C9A.2b.1 is 63 nt long and ends at 11,561 while the 1st exon of transcript Y74C9A.2b.4 is 59 nts and ends at 11,557.
So there you go: two (slightly) different transcripts.


text:
differential gene expression complex design no replicates
We have an experimental design as seen below
Where we administered drug at 0 min for each mouse genotype and took them down at given below intervals. wt and ko mouse models  were administered only with the vehicle which would be our negative control.
Performed rnaseq (total rna) truseq stranded protocol would like to see differences between each time point compared to the vehicle, used star aligner to align to mm9 , htseq for gene counts .
sample  type    time
WT  Vehicle     30 min
KO  Vehicle     30 min
WT  Drug    30min
KO  Drug    30min
WT  Drug    1 hour
KO  Drug    1 hour
WT  Drug    2 hour
KO  Drug    2hour
WT  Drug    3 hour
KO  Drug    3hour

We do not have any replicates and hence cannot use any of the differential expression statistical tools (asked this question at bioconductor where Dr. Love himself replied link to question asked) where he suggested a linear model.
we were wondering if someone could help us with any such design (papers/previously performed analysis) what we would like to see is comparison/differences (significant) between these groups
(veh vs 30min, veh vs 1 hour, beh vs 2 hours, veh vs 3 hours)

What is your design, exactly? What are we supposed to understand from that vague table? What are the times? The time you administered the drug? The length of time the drug had to act before administering an antagonist? The time you spent observing it? How is this related to bioinformatics? Please [edit] your question and explain what you are doing and what you expect from a good answer.

It's unlikely that you'll find any papers doing this because it'd be hard to get one accepted if it were largely based on such an experiment. I have to wonder why you expended time/resources on this experiment when you had to know that this was going to be a major problem.

@DevonRyan we have experimental (wet lab results supporting) this, realized it will be difficult to go ahead with rna seq only after i posted it on deseq forum (bioconductor), we are not sure if this is a dead end or we can at least see some basic differences using very basic stats

@terdon thank you for feedback i have edited the question, please let me know if further information is required

So you administered the drug and then killed the mice after the specified time. Then what? Did you perform rna seq? Of what? Whole transcriptome? What are you trying to show? What data do you have? Are you really trying to draw conclusions based on a single individual in each condition?

Although it is not recommended to use no replicates, in the edgeR manual they give some advice on how to go on with no replicates design. See page 21 of user guide. You can e.g. estimate a BCV value. Of course it is still a trick and not sound statistics.

Mike Love is right. If the response you are looking for is linear in terms of change per minute, the most productive approach is likely to be fitting a linear model. You might get something out of this because the difference between successive time points represents multiple measurements of the rate of change over time. The biggest problem is that the time points aren't independent. You would need to fit a model that estimated that rate for both WT and KO and then test the difference in rates.
I see three ways forward, in all cases you want to test the last term in the model ~time + genotype + time:genotype

Ideally you would like to fit a negative binomial linear model to the read counts. But that is going to be difficult as you have no replicates from which to assess the dispersion. Thus you could generate read counts using your favorite read count pipeline (STAR+featureCounts/Salmon/Kallisto) and take the approach in the edgeR guide noted by b.nota and then use edgeR to test the fitted model.
You could use something like salmon or kallisto to estimate the TPMs assume that log TPMs are normally distributed and model them with plain lm in R, using something like broom to model each gene separately. A tutorial by David Robinson, the author of broom on doing this sort of thing can be found here
You could use the Salmon or Kallisto bootstraps to estimate technical variation (but not biological) and then I suspect you could use Sleuth to fit the model.

In any case these results should be treated with scepticism andmostly treated as hypothesis generating. You might find genes that you want to go back and confirm and you might see some patterns of interest if you look at some enrichments etc. which you can then design experiments to test.

thank you for the information is there a tutorial for the lm and broom methods that you have mentioned would help me understand it better, as R is something I have just started using appreciate it.

I have added a link to a broom tutorial.

it is really an awesome resource many thanks!! (I would not have found this), and keep everyone posted hopefully i can pick up on using R and work this out within a couple of days.


text:
Is there a standard k-mer count file format?
I am doing a research project involving calculating k-mer frequencies and I am wondering if there is any standard file format for storing k-mer counts.

Not as far as I am aware. The Ray assembler used to (and possibly still does) store the kmers as FASTA files where the header was the count of the sequence, which I thought was a pretty neat bastardisation of the FASTA file format. It looks like this format is also used by Jellyfish when reporting kmer frequencies by the dump command (but its default output format is a custom binary format):

The dump subcommand outputs a list of all the k-mers in the file associated with their count. By default, the output is in FASTA format, where the header line contains the count of the k-mer and the sequence part is the sequence of the k-mer. This format has the advantage that the output contains the sequence of k-mers and can be directly fed into another program expecting the very common FASTA format. A more convenient column format (for human beings) is selected with the -c  switch.

Jellyfish changed their internal format between v1 and v2 (both not FASTA), because they changed to doing counts based on bloom filters. Jellyfish2 has an optional two-pass method that sets up a bloom filter intermediate file to record kmers, and multiple different final reporting formats.
Khmer also uses bloom filters, but in a slightly different way. It also has been extended to be useful for partitioning and comparing datasets.

Based on my experience reviewing k-mer counting software and as a core contributor to the khmer project, I can confidently say that there is no widely used standard format.


text:
How can I use Nanopore reads to close gaps or resolve repeats in a short-read assembly?
Low coverage MinION reads should be useful to close gaps and resolve repeats left by short-read assemblers. However, I haven't had any success with the software I know about. I'm aware of the following packages, either for scaffolding or closing gaps in short-read assemblies using long reads:

LINKS
OPERA-LG
npScarf
PBJelly
SSPACE-longread

I've tried npScarf and LINKS, and they ran successfully but didn't resolve the gaps in my assembly. I couldn't get PBJelly and OPERA-LG to run with the current versions of BLASR and samtools, and both packages seem not to be maintained. I have not tried SSPACE-longread because it's not open source.
What software can I use to fill gaps and resolve repeats in a short-read assembly with low-coverage Nanopore data?

More information
I'm finishing a mitochondrial genome. I have a 17 kb short read assembly with one gap. This was made from 2x150b paired-end reads from a TruSeq PCR-free library with a 400 b insert size. Molecular data suggests the genome may be around 32 kb.
After mapping short reads against the short read assembly, Pilon identified a tandem repeat of 156 bp in an AT-rich region, but wasn't able to close the gap.
I have 1.4 Gb of MinION rapid reads with an L50 of 8808 b. I mapped these reads against the short read assembly, and I can see reads that span the gap. I seem to have the information required to close the gap, but I don't know how to do it.
The organism has a nuclear genome > 600 Mb, so long read coverage is low. Despite this, I tried a de novo Canu assembly of the Nanopore reads, and I pulled out a 39 kb contig containing mitochondrial genes. Most genes on this contig are duplicated and fragmented, and Pilon wasn't able to improve it.
Thanks for reading!

You might want to look into Unicycler (manuscript with more information can be found here); even though it is supposed to be used with bacterial genomes only, it might work well with a small genome such as a mitochondrion's.

Beware that if you happen to have very long reads, you might end up with an assembly with multiple copies of the circular genome: you might want to look into circlator then.

Mitochondria are bacteria. Their long co-evolution as endosymbiotes may have had an influence in some of their genomic characteristics, though.

Exactly; in fact, I'm not aware of any bacteria with a ~32kb genome :)

This worked well for me. I got a single 21 kb circular contig with all the expected mt genes. Thanks for alerting me to this software!