text:
How to install DnaSP on a Mac
I need to install DnaSp, a nice GUI originally written for Windows for the population genetics of nucleotide sequence data.
However, its not working on my mac. I'm using wine to install it using these instructions for wine. Wine usually successfully ports Windows programs.
it successfully installs but then I get this error
Mollies-MacBook-Air:dnasp51001 molliepassacantando$ WINEPREFIX=~/.wine64 wine regsvr32.exe scrrun.dll mfc40.dll threed32.ocx
regsvr32: Successfully registered DLL 'scrrun.dll'
regsvr32: Successfully registered DLL 'mfc40.dll'
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
err:typelib:sltg_get_typelib_ref Unable to find reference
regsvr32: Successfully registered DLL 'threed32.ocx'
Mollies-MacBook-Air:dnasp51001 molliepassacantando$ winedevice.exe(30828,0x405ed000) malloc: *** error for object 0x401cea08: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug

and cannot run it.
Has anyone successfully installed DnaSp on a Mac or know of an alternate program I could use? I need to be able to calculate haplotype frequencies.

Hey, I see you found your way here. Great! But as I asked you on [ubuntu.se], why are you using wine? The DnaSP authors state that on a mac  "You can run DnaSP using SoftWindows, Virtual PC, Parallels Desktop or VMware Fusion software [sic]". Also, please [edit] your question and explain the kind of analysis you need to do. The vast majority of bioinformatics tools don't work on Windows and work fine on a UNIX-based system like Mac. There are almost certainly other choices but we need to know what you are trying to achieve.

Welcome to Bioiformatics stack exchange. We can not suggest you an alternative if we do not understand your problem. Could you [edit] your question to elaborate bit more what you want to use the software for?

@kamil i need it to calculate haplotype frequencies, im using wine because i wasnt sure what other platform to use, and i dont want to pay for software, as i found you needed to do with the suggested programs

@MolliePassacantando that sounds reasonable, I bet there will be a lot of free softwares doing that. However, please [edit] your question so people do not have to search this information in comments. Also do not hesitate to add more information about software you already tried...

@terdon i have edited the question, and am not sure why i cannot run the program with a mac. the file i download from the website is a .msi, which i open using wine, because mac does not read that format

This looks like a WINE error unrelated to DnaSp. How did you install it?

youre right! it was a WINE error, I managed to get the error to stop showing up, but im not sure how. but it works!

From the comments @MollePasscantando successfully ported DNAsp onto a Mac (Air) via Wine, but was not sure how. This is good news because Wine is very commonly used for porting Windows exes. The website states,

Oracle VirtualBox (free open-source). Parallels, CrossOver, VMware
Fusion (paid-versions; trial period)

From the website it was update in February 26, 2019 and its current Beta Version is 6.12.03.
There was a major published update to the package in Mol Biol Evol in 2017 (version 6) for large data sets (NGS).
Thus the question is still relevant because they still have not written an OSX specific version and this will limit its performance for the data sets it is now targeted for.


text:
Reordering scaffolds according to a reference without a genetic map
I am trying to reorder scaffolds of a rice species, but no genetic map is available right now. Oryza sativa Japonica is a close relative of this rice species. Mummer was used to do a whole genome alignment, and I am trying to reorder scaffolds according the mummer result. But I found it's too complicated to do this by myself, many scaffolds may align to multiple Oryza sativa Japonica chromosomes or multiple loci in the same chromosome.
What software can I use to handle this task?

I've been using ATAC to successfully align different versions of the rice genome assembly. There is a chance it will work for you between two closely related species (potato and tomato work well for example).
It's not massively user friendly, but it's very fast, and it produces long, 'clean' alignments (i.e. almost all 1 to 1 alignments). The output format is a bit cryptic (zero-based start position and length instead of the more usual 1-based start and end positions), but as I said, it should basically output the order for you in about 30 minutes.

That looks like a good solutions. I hopt this will help me to align scafolds

Have you tried Mauve alignment? Its pretty easy once you become familiar and has a GUI for further ease of use. Additionally there are a few online tutorials on how to re-order contigs/ scaffolds using this software. Heres one I use when in need.
Mauve contains a function called Mauve Contig Mover (MCM) which can be used to a) compare an assembly to a reference and/or b) rearrange contigs to improve assembly quality.
I use the tutorial as above which uses the GUI for Mauve. I

Could you please expand this and give an example of how Mauve can be used for the OP's purposes?

source: Stack Exchange

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

sort: votes


text:
tools to reconcile experimental transcripts with reference annotation
Looking for tools to reconcile alignment file of experimental transcripts mapped to genome (SAM/BAM) with the reference transcriptome annotation (GTF) from Ensembl (organism: D. melanogaster).
The aim would be to check which transcripts reported in the alignment file are also reported in the reference annotation, generating a new annotation (e.g., GTF) including the new transcripts as seen in the alignment file (experimental transcripts to genome).
Any idea?
-- Edit:
Up to now, I tried an in-house script that does the following:

generate a BED file from the BAM file with the exon coordinates (start/end) per transcript - using BEDOPS v2p4p27
convert the reference GTF to BED to retrieve the coordinates per annotated feature
compare the experimental vs the reference BED intervals
infer the suite of annotated exons per experimental transcript, thus, ideally, the isoform it corresponds to

However, with this simple approach, I only get around 20% of the experimental transcripts that can be identified in all their length using reference annotation. This might be due either to a too simple approach (e.g. if the start/end of the exon in the experimental transcript is shifted by 1 base as compared to the reference, it won't be detected), or to a real biological signal (since I am also using data from long-read sequencing technologies).

Can you convert the GTF file (from the BAM file via BED) to GTF? Then you can use cuffcompare. I've never tried making a GTF out of BAM files, so I haven't a clue if that works.

I will try that, and post an answer here if I get to a solution

I've never tried this myself, so I don't know how easy this is...
One option would be to start with GMAP, which is meant to align whole transcripts against the genome. The really nice thing about this is that it can directly produce GFF3 files. You can then use that with your Ensembl GTF with cuffcompare or whatever the equivalent is in stringTie. You should then be able to get the information you want from the transfrag codes.

I think that's exactly the solution we are converging to ;)

As per my answer to @_julien_roux on twitter:
Trying to find novel transcripts within the context of an existing annotation is much less straightforward. You probably need to do a "genome-guided assembly" with Trinity and PASA:
http://pasapipeline.github.io/#A_ComprehensiveTranscriptome
We did something similar in much simpler organisms in our recent paper:
http://dx.doi.org/10.1186/s12864-017-3505-0
Your biggest problems are going to be alternative spliced gene models, the large amount of fragmentation inherent in de novo transcriptome assembly and largely underreported non-coding RNAs in most genomes. It was easier for us as the dictyostelids have few and small introns - fly is going to be a different ball-game. Be prepared to spend a lot of time optimising parameters before finding a result you're happy with (maybe) ;)
You may think the fly genome pretty good and complete, but be prepared to find lots of errors, corrections and novelty. We were very surprised at the amount of new annotation we were able to find in D. discoideum - an organism where every single gene had been manually curated!
Good luck!

Sure, we know this is not a simple project ;) We now have the reconstructed transcripts, from different methods and types of data and would like to evaluate them by comparing to reference annotation (one out many ways we are trying to evaluate these)

@jul635 ah ok. Then I'd strongly recommend PASA, it gives an interactive report of all the changes made to an existing annotation by your transcripts. I'd also recommend Transrate as great tool for quality assessing your assembled transcripts.

Thanks to @Devon_Ryan's suggestion, I had a go at using GMAP with our own de-novo assembly. Here's what I started with:

Genome assembled from nanopore long reads, corrected using Illumina cDNA reads
genome-guided Trinity transcriptome assembly, with reads grouped based on mapping to the assembled genome, then assembled by Trinity

The transcriptome assembly was filtered to only include transcripts that had reasonable coverage from cDNA reads, and further filtered to cluster similar transcripts and only use the longest isoform. This process is detailed in our paper. If I were to do that again, I'd probably create a supertranscript containing all exons, rather than the longest isoform.
I first created an index for the genome using gmap_build:
gmap_build -k 15 -j 8 -d 'MINBRA' \
  Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta -D MINBRA_GMAP

I then used gmap to map the transcriptome to the genome, generating gene models as a GFF3 file:
gmap -t 10 -d MINBRA -D ./MINBRA_GMAP \
  -A locusTagged_CDLI_CD98LMOHC50_TBNOCFED.tran.fasta \
  -f 2 > Nb_TrinTran_vs_Canu.gff3

Here's an example from the annotation file that was produced by GMAP:
...
tig00022147 MINBRA  gene    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path1;Name=MINBRA_TGG_483_c1_g1_i1
tig00022147 MINBRA  mRNA    413327  414778  .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path1;coverage=100.0;identity=99.9;matches=1452;mismatches=0;indels=2;unknowns=0
tig00022147 MINBRA  exon    413327  414778  99  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00022147 MINBRA  CDS 413417  413557  100 +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna1.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna1;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00001947 MINBRA  gene    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.path2;Name=MINBRA_TGG_483_c1_g1_i1
tig00001947 MINBRA  mRNA    37526   38961   .   +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.path2;coverage=100.0;identity=98.7;matches=1435;mismatches=1;indels=18;unknowns=0
tig00001947 MINBRA  exon    37526   38961   98  +   .   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.exon1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 1454 1 .
tig00001947 MINBRA  CDS 37614   37750   97  +   0   ID=MINBRA_TGG_483_c1_g1_i1.mrna2.cds1;Name=MINBRA_TGG_483_c1_g1_i1;Parent=MINBRA_TGG_483_c1_g1_i1.mrna2;Target=MINBRA_TGG_483_c1_g1_i1 231 91 .
###
tig00022147 MINBRA  gene    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.path1;Name=MINBRA_TGG_489_c13_g1_i2
tig00022147 MINBRA  mRNA    464550  465456  .   +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.path1;coverage=100.0;identity=99.3;matches=907;mismatches=0;indels=6;unknowns=0
tig00022147 MINBRA  exon    464550  465456  99  +   .   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.exon1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 913 1 .
tig00022147 MINBRA  CDS 464909  465053  98  +   0   ID=MINBRA_TGG_489_c13_g1_i2.mrna1.cds1;Name=MINBRA_TGG_489_c13_g1_i2;Parent=MINBRA_TGG_489_c13_g1_i2.mrna1;Target=MINBRA_TGG_489_c13_g1_i2 507 361 .
###
...

GMAP told me on standard output what transcripts were not mappable to the genome (136 of 36830, or about 0.4%).
I've now used STAR to map these cDNA reads to the transcripts (using the GFF3 annotation file for initial evidence) in order to confirm that the mapping carried out by GMAP is reasonable. I ended up throwing out another 7.4% of transcripts. First, I generated a genome index using the generated gff3 file:
~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
   --runMode genomeGenerate --genomeDir ./star_genome \
   --genomeFastaFiles Nb_ONTCFED_65bpTrim_t1.contigs.hpcleaned.fasta \
   --sjdbGTFfile Nb_TrinTran_vs_Canu.gff3 \
   --sjdbGTFtagExonParentTranscript Parent

Followed by mapping my RNASeq experiments using STAR. STAR can sort the genome mapping, but not the transcriptome mapping, so an additional sort/index was done on that:
for x in $(seq 1 9); do
  ~/install/star/STAR-2.5.4b/bin/Linux_x86_64/STAR --runThreadN 10 \
     --readFilesCommand zcat --genomeDir ./star_genome \
     --readFilesIn *-1563-0${x}-*.fastq.gz \
     --outFileNamePrefix mapped/1563-0${x}_vs_BNOCFED \
     --outSAMtype BAM SortedByCoordinate --quantMode TranscriptomeSAM
  samtools sort mapped/1563-0${x}.toTranscriptome.*.bam \
    > mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
  samtools index mapped/1563-0${x}.toTranscriptome.sortedByCoord.out.bam
done

BAM hits were collated into a single file:
(for x in $(seq 1 9);
    do samtools idxstats 1563-0${x}.toTranscriptome.sortedByCoord.out.bam | \
      perl -pe "s/^/1563-0${x}\t/";
 done) | gzip  > readCounts_all_vs_transcriptome.tsv.gz

Which was then processed using R to identify unmapped transcript annotations:
#!/usr/bin/Rscript
data.df <-
    read.delim("readCounts_all_vs_transcriptome.tsv.gz", header=FALSE,
               col.names=c("Sample","Transcript",
                           "Length","Mapped","Unmapped"));
gene.counts <- tapply(data.df$Mapped + data.df$Unmapped,
                      list(data.df$Transcript, data.df$Sample), sum);
sum(rowSums(gene.counts)==0) / nrow(gene.counts);
TranscriptsToRemove <- rownames(gene.counts)[rowSums(gene.counts) == 0];
cat(TranscriptsToRemove, file="transcripts_to_remove.txt", sep="\n");

The unmapped annotations were then excluded from the gff3 file:
grep -v -F -f mapped/transcripts_to_remove.txt Nb_TrinTran_vs_Canu.gff3 \
  > mappable_Nb_TrinTran_vs_Canu.gff3

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


text:
How can I do an overlapping sequence count in Biopython?
Biopython's .count() methods, like Python's str.count(), perform a non-overlapping count, how can I do an overlapping one?
For example, these code snippets return 2, but I want the answer 3:
>>> from Bio.Seq import Seq
>>> Seq('AAAA').count('AA')
2
>>> 'AAAA'.count('AA')
2

For Biopython 1.70, there is a new Seq.count_overlap() method, which includes optional start and end arguments:
>>> from Bio.Seq import Seq
>>> Seq('AAAA').count_overlap('AA')
3
>>> Seq('AAAA').count_overlap('AA', 1, 4)
2

This method is also implemented for the MutableSeq and UnknownSeq classes:
>>> from Bio.Seq import MutableSeq, UnknownSeq
>>> MutableSeq('AAAA').count_overlap('AA')
3
>>> UnknownSeq(4, character='A').count_overlap('AA')
3

Disclaimer: I co-contributed the .count_overlap() methods with Peter Cock, see 97709cc

I've encountered this problem before, and used python re module to solve this problem.
import re
all = re.findall(r'(?=(AA))','AAAA')
counts = len(all)

You can get more details in this thread

This is a valid approach, and you could add start and end arguments easily enough, but it's not very memory efficient for large sequences. re.finditer would be better for that. Also, UnknownSeq type objects require a rather different approach to maximize memory efficiency

You can use finditer from python's re module. The advantage of this approach is it allows for getting the indices of those matches, which could be handy down the track.
>>> import re
>>> matches = re.finditer(r'(?=(AA))', 'AAAA')
>>> indices = [match.span(1) for match in matches]
>>> indices
[(0, 2), (1, 3), (2, 4)]
>>> num_matches = len(indices)
>>> num_matches
3


text:
Minfi returning incorrect beta values
UPDATE: I found the solution. I was using normalized values and GEO was using raw beta values.
I'm trying to link GEOquery and minfi.
Specifically I want to obtain beta values from the idat files on GEOquery. I was following this guide: https://kasperdanielhansen.github.io/genbioconductor/html/minfi.html, up until the preprocessing part. Meaning that I was able to obtain the RGset. Then I used my own preprocessing code to obtain the beta values. However, I cross checked them with the beta values that were on GEO and they weren't consistent.
For example, the accession number I used was GSE68777. So I went to that study on GEO and clicked on the first sample: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1681154. Then I scrolled down and clicked "Download full table" to download the samples and beta values in a text file.
Then I typed head(beta) and chose the first sample. Then I did command F in the text file for that sample and it's value there wasn't the same as the value from the beta data table. Hopefully you can help find the error.
Here is the code I'm using:
library(GEOquery)
library(minfi)
library("IlluminaHumanMethylation450kanno.ilmn12.hg19")
library("IlluminaHumanMethylation450kmanifest")

######## Code copied from 450k Guide ########

getGEOSuppFiles("GSE68777")
untar("GSE68777/GSE68777_RAW.tar", exdir = "GSE68777/idat")
head(list.files("GSE68777/idat", pattern = "idat"))
idatFiles <- list.files("GSE68777/idat", pattern = "idat.gz$", full = TRUE)
sapply(idatFiles, gunzip, overwrite = TRUE)
rgSet <- read.metharray.exp("GSE68777/idat")
geoMat <- getGEO("GSE68777")
pD.all <- pData(geoMat[[1]])
pD <- pD.all[, c("title", "geo_accession", "characteristics_ch1.1", "characteristics_ch1.2")]
names(pD)[c(3,4)] <- c("group", "sex")
pD$group <- sub("^diagnosis: ", "", pD$group)
pD$sex <- sub("^Sex: ", "", pD$sex)
sampleNames(rgSet) <- sub(".*_5", "5", sampleNames(rgSet))
rownames(pD) <- pD$title
pD <- pD[sampleNames(rgSet),]
pData(rgSet) <- pD

Could you actually post the numeric values, so I don't have to rerun and do all again? Did you check if for the same position do you have the same beta value? Maybe you weren't comparing the same positions

Presumably the authors used a specific preprocess*() function. Since they didn't put it in the paper, you'll need to ask them.

You really need to show us the exact value you’re getting, and where you think you’re getting it from (which probe). This sounds like it’s just normal floating point imprecision. Meaning the value is present in the text file but with different rounding.

I figured it out. They used the raw beta values but I was using normalized values. Thanks for the comments though

@EzraBekele Could you write a short answer stating this? It makes it more easily findable for other users.


text:
Profile of conservation between bacterial sequences and human protein
My goal is to make a conservation plot between bacterial sequences and a human protein.
So far, I have a FASTA file of the protein, and a FASTA file with the sequences of the proteins from the BLAST results.
In my initial attempts, I have tried to do this using msa software:
I have downloaded clustalx and run a profile alignment using the single protein sequence as Profile 1 and the sequences from the BLAST results as Profile 2, and selected the option "Align Sequences to Profile 1".
I'm having difficulty interpreting and dealing with the results. Since the majority of the sequences aren't aligned, the results are a total mess, with a huge number of dashes added.
My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.
ie:
I would like to have a graph with my protein on the x-axis and a plot like the regions of high conservation for the multiple alignments, except with the spikes corresponding to a high number of BLAST hits. This would allow me to identify regions of my protein that have higher sequence similarity to bacteria than others.
Is there a better way to achieve this or to salvage the results from the profile alignment?
Thanks.

I curious why you want to use a profile alignment. Wouldn't it be easier to just add your sequence to the sequences from the blast results and run a standard multiple sequence alignment?

Also, could you clarify what you're trying to achieve? You say you want the number of blast hits per amino acid, but the figure you show is plotting amino acid conservation. Do you care if the amino acids are the same at a given position or just if the high-scoring segment pairs overlap the position?

We could help much more if you gave us the alignment in text form or even a few of your sequences so we can align them here. The image is informative but we can't really work with it. That said, please do explain what you're trying to do here. I am having trouble understanding what the number of blast hits per amino acid would convey. Why would you need that value? And are you sure you want all of those sequences? If the majority are not aligned, are you sure they're relevant? What is the actual question you are trying to answer using this alignment?

A profile alignment generates a new alignment from two sets of aligned sequences (the "families"), where one particular case is between one family and one sequence. Therefore to run clustalx, you need to first align the proteins from BLAST results, before using them as your Profile 2. The other questions remain, about number of blast hits or final goal...maybe there is some confusion with the position-specific score matrix from PSI-BLAST, which also generates a score?

Is there a way to visualize the position specific scoring matrix generated from a PSI-BLAST with a line plot that would resemble the graph from the msa in my original post? My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

@bluescholar1212 have you looked at Jalview? It is designed as an MSA viewer, plus has many ways for viewing conservation which you can extend yourself.

@ithinkiam I have taken a look at Jalview and it looks like it produces some great visualizations. Would a BLAST result be a valid Jalview input?

AFAIK, no it wouldn't be a valid input, but you could use the blast ids to make the alignment and provide the multiple alignment to Jalview. Also please consider [edit]ing the question, you want to have a profile of conservation between bacterial sequences and the human protein, isn't it?

@Llopis thanks, I went ahead and edited my original question, and am going to add a better description as well shortly.
Regarding how best to achieve this, am I better off using a different approach, or is that a reasonable workflow to blast, align the results, and visualize with Jalview? I think potentially a line plot showing consensus percentages of the results against the human protein would be an appropriate metric.

I'm not sure what you mean by it being a mess? That looks like a pretty good alignment to me, and it most certainly has aligned all the sequences. You are nearly always going to have gaps (see my MSA below, which are all paralogs from the same genome so they couldn't really be more closely related, and yet, there are gaps). That comes with its own set of problems mind you, as you need to decide how you're going to deal with them.
You also don't necessarily have to do a profile alignment. Those sequences don't look massively divergent to me, so straight forward sequence alignment would probably give you more or less the same result.

My goal is to visualize the number of BLAST hits per amino acid in my protein of interest.

This doesn't really make sense. You should really be searching for the number of BLAST hits with a particular domain, if you want to infer homology.
If you want a per-position visualisation of the conservation, you could look at the shannon entropy for each column and plot that. I wrote a script to do just that a little while ago: https://github.com/jrjhealey/bioinfo-tools/blob/master/Shannon.py
Just beware it's not super well tested yet. Feed an MSA in with as many sequences as you want to analyse, but you'll have to have identified the sequences and done the alignment first.
For example, given this MSA:
    16    149
PAU_02775  MSTTPEQIAV EYPIPTYRFV VSLGDEQIPF NSVSGLDISH DVIEYKDGTG
PLT_01696  MSTTPEQIAV EYPIPTYRFV VSIGDEQIPF NSVSGLDISH DVIEYKDGTG
PAK_02606  MSTTPEQIAV EYPIPTYRFV VSIGDEQVPF NSVSGLDISH DVIEYKDGTG
PLT_01736  MSTTPEQIAV EYPIPTYRFV VSIGDEKVPF NSVSGLDISH DVIEYKDGTG
PAK_01896  MTTTT----V DYPIPAYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PAU_02074  MATTT----V DYPIPAYRFV VSVGDEQIPF NSVSGLDITY DVIEYKDGTG
PLT_02424  MSVTTEQIAV DYPIPTYRFV VSVGDEQIPF NNVSGLDITY DVIEYKDGTG
PLT_01716  MTITPEQIAV DYPIPAYRFV VSVGDEKIPF NNVSGLDVHY DVIEYKDGTG
PLT_01758  MAITPEQIAV EYPIPTYRFV VSVGDEQIPF NNVSGLDVHY DVIEYKDGIG
PAK_03203  MSTSTSQIAV EYPIPVYRFI VSIGDDQIPF NSVSGLDINY DTIEYRDGVG
PAU_03392  MSTSTSQIAV EYPIPVYRFI VSVGDEKIPF NSVSGLDISY DTIEYRDGVG
PAK_02014  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAU_02206  MSITQEQIAA EYPIPSYRFM VSIGDVQVPF NSVSGLDRKY EVIEYKDGIG
PAK_01787  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PAU_01961  MSTTADQIAV QYPIPTYRFV VTIGDEQMCF QSVSGLDISY DTIEYRDGVG
PLT_02568  MSTTVDQIAV QYPIPTYRFV VTVGDEQMSF QSVSGLDISY DTIEYRDGIG

           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINISLRKGV FSGDTKLFDW INSIQLNQVE KKDISISLTN
           NYYKMPGQRQ AINITLRKGV FSGDTKLFDW LNSIQLNQVE KKDISISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NHYKMPGQRQ LINITLRKGV FPGDTKLFDW LNSIQLNQVE KKDVSISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NYYKMPGQRQ SINITLRKGV FPGDTKLFDW INSIQLNQVE KKDIAISLTN
           NWFKMPGQSQ LVNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NWFKMPGQSQ STNITLRKGV FPGKTELFDW INSIQLNQVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NYYKMPGQIQ RVDITLRKGI FSGKNDLFNW INSIELNRVE KKDITISLTN
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPTITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD
           NWLQMPGQRQ RPSITLKRGI FKGQSKLYDW INSISLNQIE KKDISISLTD

           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL ISPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TSPSFDATSN EVAVQEITLT ADRVTIQAA
           EAGTEILMTW SVANAFPTSL TAPAFDATSN EVAVQEISLT ADRVTIQAA
           ETGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EVGTEILMTW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVTIQAA
           EAGTEILMSW SVANAFPTSL TSPSFDATSN DIAVQEIKLT ADRVMIQAA
           ETGSQILMTW NVANAFPTSF TSPSFDAASN DIAIQEIALV ADRVTIQAP
           EAGTEILMTW NVANAFPTSF TSPSFDATSN EIAVQEIALT ADRVTIQAA
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLT ADRVIMQAV
           DAGTELLMTW NVSNAFPTSL TSPSFDATSN DIAVQEITLM ADRVIMQAV
           DTGSEVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           DTGSKVLMSW VVSNAFPSSL TAPSFDASSN EIAVQEISLV ADRVTIQVP
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEMSLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH
           ETGSNLLITW NIANAFPEKL TAPSFDATSN EVAVQEISLK ADRVTVEFH

You'd get this plot:

I would like to have a graph with my protein on the x-axis and a plot
  like the regions of high conservation for the multiple alignments,
  except with the spikes corresponding to a high number of BLAST hits.
  This would allow me to identify regions of my protein that have higher
  sequence similarity to bacteria than others.

I'm not sure your logic is quite right here though. Blast won't give you hits depending on a particular position. It's a local aligner, so it'll just return you hits where at least some part of your query matches at least some part of another.
What you could do is take the logic in the script above, and just use a different metric. For example, perhaps you could count the proportion of sequences which have the most common amino acid at a given position within your MSA. That would be fairly crude though.
As you say in the comments,

My final goal is to produce a plot visualizing regions of high bacterial sequence similarity to my human protein of interest.

your original alignment will show you this intrinsically, if only you include the sequences of all the BLAST hits in the first place. Thus your work flow will be:

Blast sequence of interest.
Download all/as many hits as you want (bear in mind the E-value/bitscore and the number of hits you get. It might only be a few dozen, in which case you can use the lot, but if not, just take all the hits below a certain cut-off.)
Align all the sequences.
Look at the column scores for the whole MSA. You can use whatever metric of conservation you like really. Might be as simple as proportion of sequences with the most common residue, or something more complex like Shannon entropy (though as you can see in the graph above Shannon entropy can be kinda noisy) etc.

this is an excellent answer. I think that this is the work flow that I will adopt. While researching this I also found this blog post link with a script that produces a pretty smooth conservation plot using the Savitzky-Golay technique. The script seems to follow a similar workflow to what you proposed. Now I will need to work on which blast parameters and metrics of conservation will be appropriate for this project. Thanks for your help.

If you don't mind abstracting the data slightly, the script can calculate a moving average/running mean which smooths the data a bit.

this may be worth posting a new question to answer, but in the link that I posted, the author mentions you could easily plug in a new function for the remote_blast function depending on usage. I have a list of parameters I would like to use for my blast, but ideally I would like to access the remote database. Do you know how I could change the function in that script to specify more specific parameters that I need?

I'm not sure what you mean exactly? With the python implementations of remote blast you should have most, if not all, the options of a standard blast. What sort of parameters are we talking about? It's easy enough to add additional arguments to the argparse framework if you want to specify blast thresholds etc. It's pretty much just a case of adding in another function or two to do the heavy lifting as far as generating the MSA for you goes, and then that can be passed directly to the other functions of the file.


text:
How to interpret methylation calls from Bismark on opposite strands?
I'm looking at Reduced representation bisulfite sequencing (RRBS) data from ENCODE, and to align the FASTQ files I've used Bismark with Bowtie 1. When I load the resulting BAM file into R with Rsamtools and GenomicRanges, I get something like this:
> reads
GRanges object with 2 ranges and 5 metadata columns:
      seqnames         ranges strand |                     seq                                   XM          XR          XG        NM
         <Rle>      <IRanges>  <Rle> |          <DNAStringSet>                          <character> <character> <character> <integer>
  [1]     chrM [16523, 16558]      - | ATAAAACCTA...CCCTTAAATA .....h........h.......z.............          CT          GA         3
  [2]     chrM [16524, 16559]      + | TAAAGTTTAA...TTTTAAATAA .....hh.......hhh.h.z...hhhh........          CT          CT        11
  -------
  seqinfo: 25 sequences from an unspecified genome

For the read on the - strand (which is aligned to the "G->A"-converted reference), how should I compare methylation calls to those on the + strand, since they don't line up?
For example, the z in the methylation string of the - strand is one position ahead of the z in the + strand (which makes sense because of symmetric CpG methylation). But how should I determine whether these two methylation calls are "essentially the same" or not?

The strand the bismark reports is related to the strand from which the read originated, not necessarily how it's aligned. So, alignments on the + strand shouldn't have calls overlapping those on the - strand, since you can't have a C in the same place on the same strand. One should often see Z/z next to each other on opposite strand, like in your example, since these are CpG (so it should be a Z/z on the + strand and then a Z/z on the following base on the - strand). Thus, in your example you have two reads supporting unmethylation for the CpG as a whole (1 read for each of the Cs).
The most confusing thing about BSseq is that reads (or paired-end reads) are only ever informative for a single strand. This is actually true for all sequencing with Illumina instruments (single-stranded fragments are loaded, after all), but since strands are almost always complementary we can usually ignore this fact.
As an aside, you might find MethylDackel useful (full disclosure, I wrote it). It'll be much faster at extracting methylation calls than bismark and supports nice things like excluding regions of bias methylation and likely variant positions.

OK, that's what I wasn't sure about, because of the description given by the Babraham Institute about how Bismark works (pg. 5, for anyone who's curious). So not having the Z/z's overlap between reads on opposite strands makes sense.
So to say that "these reads are concordant in their methylation calls", would I just shift the positions that I compare between + and - strand reads by 1 to properly align them? Only for CpG methylation, not CHH- or CHG-context methylation

Yes, you can do that for CpGs.


text:
Filter Trinity transcriptome based on RNASeq reads
I have recently generated a genome-guided transcriptome with Trinity, and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is with the goal of trying to reduce the initial dataset down to something a bit more manageable (I have about 300k transcripts covering 250Mb in total, but would prefer about 1/10 of that number).
I see that Trinity has a workflow using Bowtie2 for evaluating the quality of assembled transcripts, but Bowtie2 is getting a bit old in the mapping world now. Are any of the super-fast transcript mappers (e.g. Salmon or Kallisto) appropriate to use for transcript filtering?

It's encouraged for Stack Exchange. There's even a FAQ Answer about it. It serves as a documentation of a tricky problem/solution I have, and also allows others to contribute and offer other solutions that I wasn't aware of (like your own).

You may try using hisat2 instead of bowtie2 in the workflow you mention: hisat2 almost has the same user interface as bowtie2, and could be faster.

I've tried both HISAT2 and Bowtie2 for mapping. HISAT2 is splice-aware, so should be better for mapping RNASeq reads to a perfect genome, but Bowtie2 gives slightly better BUSCO scores for sequencing error correction based on mapping (in my particular situation).

I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:

Run BUSCO in short/transcript mode on the Trinity-generated sequences
Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
Sub-select the transcripts based on this threshold, using the Salmon results table

Step 1 -- BUSCO
python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10

Step 2 -- Salmon
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A

Step 3 -- Results merge (using R)
#!/usr/bin/Rscript
busco.table <-
    read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
               header=FALSE, stringsAsFactors = FALSE,
               col.names=c("Busco id","Status","Sequence","Score","Length"));
busco.table <- subset(busco.table, Status != "Missing");

count.table <- read.delim("quant/reads-all_quant_ISR/quant.sf",
                          stringsAsFactors=FALSE);

busco.expr.df <- merge(busco.table, count.table,
                       by.x="Sequence", by.y="Name", all.x=TRUE);

Step 4 -- Count thresholding (visual exploration)
pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
options(scipen=10);
busco.numReads <- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                              function(x){
                                  exp(mean(log(x[x>0])))}), decreasing=TRUE);
par(mfrow=c(1,2));
plot(x=seq(0,1,length.out=length(busco.numReads)),
     busco.numReads, log="y", ylab="Number of Mapped Reads",
     yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
     xleft=0.9, xright=par("usr")[2],
     col="#80808020", border=NA);
plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
     busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
     xlab="Proportion of Complete BUSCO sequences");
axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
invisible(dev.off());

Looks like a value of 50 will be suitable.
Step 5 -- Sub-select the transcripts
cat(file="HighCount_Transcripts_Num50.txt",
    subset(count.table, NumReads >= 50)$Name, sep="\n");

Then filter using a suitable FASTA selection utility:
pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta

I really like the apporach! How do you handle/treat replicates?

Do you mean identical transcripts (or sub-transcripts)? I've used a very strict cd-hit-est clustering to clean up a few (98% identity), but the BUSCO results still indicate that there's a lot of duplication going on. I'm hesitant to use a looser identity threshold because I don't want to include any of the real gene duplications.

I would recommend Transrate. It generates a score for each transcript based on many metrics (read their paper here) and, using that score, gives a predicted cut-off to produce the optimal transcript set.

DETONATE offers a pipeline for filtering in this way as well. But the aligner used by default is Bowtie2...

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


text:
how is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation?
For DNA/RNA quantification machines like the Bioanalyzer or TapeStation, the DNA Integrity Number (DIN) or RNA Integrity Number (RIN) numbers are quoted as a measure of the fragmentation of the material.
How is the DNA Integrity Number (DIN) calculated in Bioanalyzer/TapeStation? Is there any open-source software to parse/collate the data from the machines?

Note that DIN (DNA Integrity Number) and RIN (RNA Integrity Number) are different scores. I've had trouble in the past finding code or formulas for calculating both of these scores, which is very frustrating considering how frequently they are used in research papers and for NGS QC. The closest that I have been able to get to this are demonstrative graphs for various RIN numbers and DIN numbers, with a description of the variables that are included in the model.
The code that @719016 provided appears to attempt to calculate the RIN, not the DIN.
Just in case it's helpful, here is an extract of that code, only including lines that seem to be directly involved in calculating the actual B_RIN values. This is not runnable, because it depends on particular input files, and it's still pretty unintelligible:
myfiles<-list.files()
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str3<-paste0(substr(targets2[1,2],1,67),"_S")
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
    ##mod chip numb
    targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
    str3<-paste0(substr(targets2[1,2],1,67),"_S")
    samp_num_chip<-length(grep(str3,myfiles))
    ##mod sample numb
    str1<-substr(targets2[1,2],1,67)
    for(j in 1:samp_num_chip)
    {
        iii<-j
        str5<-paste0(str1,"_Sample",iii,".csv")
        samples1[j,i]<-str5
    }
}
samples1<-samples1[grep("2100",samples1)]
##RIN Calculator#######refined peakfinder####################
##read bioanalyzer data into a matrix called dta
##since the total RNA and mRNA assays run differently
##skip more lines for the total RNA assay
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
    x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
    dta[,i] <- x[,2]
}
time<- x[,1]
for(samples in 1:length(samples1))
{
    ##store max fluor in max.peaks matrix
    lad<-samples
    ti1<-17
    ti2<-65
    p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
    peaks <- matrix(NA, nrow=p.row, ncol=1)
    max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
    for (i in 1:length(seq(ti1,ti2,.05)))
    {
        ii<-seq(ti1,ti2,.05)[i]
        x0<- ii
        x0<-round(x0,digits=3)
        y0<- ii
        x<- ii+.05
        x<-round(x,digits=3)
        y<- ii
        max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
    }
    ##human or mouse
    cent<-max(max.peaks)*.6
    ##flatworm
    ##cent<-max(max.peaks)*.05
    ##define the middle of the segments
    max.avg<-seq(ti1+.025,ti2+.025,.05)
    ##find time associated with peaks > 5 fu and fill seconds matrix with them
    ##12.5% of max fu incase chip is fu
    sizes <- which(max.peaks>cent)
    gt35<-which(max.avg[sizes]>35)
    lt55<-which(max.avg[sizes]<55)
    tfr<-which(duplicated(c(gt35,lt55))==TRUE)
    area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],
                max.peaks[sizes][c(gt35,lt55)[tfr]],
                max.peaks[sizes+1][c(gt35,lt55)[tfr]])
    time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],
                 max.avg[sizes][c(gt35,lt55)[tfr]],
                 max.avg[sizes+1][c(gt35,lt55)[tfr]])
    t1<-time2[which(area[,1]<cent),1][1]
    t2<-time2[which(area[,1]<cent),1][2]
    t3<-time2[which(area[,3]<cent),3][1]
    t4<-time2[which(area[,3]<cent),3][2]
    s18a<-ifelse(is.na(t1,-30,round(t1,digits=1)))
    s18b<-ifelse(is.na(t3,-35.05,round(t3,digits=1)))
    s28a<-ifelse(is.na(t2,-30,round(t2,digits=1)))
    s28b<-ifelse(is.na(t4,-30.05,round(t4,digits=1)))
    eta<-which(time==s18a)
    etb<-which(time==s18b)
    tea<-which(time==s28a)
    teb<-which(time==s28b)
    s18<-sum(dta[eta:etb,lad])
    s28<-sum(dta[tea:teb,lad])
    qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
}

Mostly based on that last line, the RIN seems to be associated with the ratio of the 18s peak to the 28s peak.

The best I could find by looking at github.com is:
https://github.com/brianfleharty/Agilent-BioAnalyzer-RIN-R-code-for-Fly/blob/master/RNA_RIN.R

#set the working directory
dir<- "H:\\B_RIN\\emd"
setwd(dir)
myfiles<-list.files()
which(substr(myfiles,69,79)=="Results.csv")
results<-myfiles[which(substr(myfiles,69,79)=="Results.csv")]
n <- which(substr(myfiles,69,79)=="Results.csv")
targets <- read.delim(myfiles[n[1]], sep=",",header=F,skip=14,as.is=T)
length(n)
targets2<- read.delim(myfiles[n[1]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str1
str2<-"_S"
str3<-paste(str1,str2,sep="")
str3
chip_num<-length(n)
samp_num_chip<-length(grep(str3,myfiles))
samp_num<-length(which(substr(myfiles,69,74)=="Sample"))
samples1<- matrix(NA, ncol=chip_num, nrow=samp_num)
for(i in 1:chip_num)
{
#mod chip numb
targets2<- read.delim(myfiles[n[i]], sep=",",header=F,skip=0,as.is=T)
str1<-substr(targets2[1,2],1,67)
str2<-"_S"
str3<-paste(str1,str2,sep="")
samp_num_chip<-length(grep(str3,myfiles))
#mod sample numb
str1<-substr(targets2[1,2],1,67)
str4<-"_Sample"
for(j in 1:samp_num_chip)
{
iii<-j
end.st<-".csv"
str5<-paste(str1,str4,iii,end.st,sep="")
samples1[j,i]<-str5
}
}
samples1
grep("2100",samples1)
samples1<-samples1[grep("2100",samples1)]
samples1
write.table(samples1,file=paste(dir,"\\ba_lane.txt",sep=""), sep="\t",col.names=NA)
#RIN Calculator#######refined peakfinder####################
#read bioanalyzer data into a matrix called dta
#since the total RNA and mRNA assays run differently
#skip more lines for the total RNA assay
ramp <- colorRamp(c("darkmagenta","white"))
elec <- read.delim("ba_lane.txt", sep="\t",header=T,as.is=T)
qc.mat <- matrix(NA, ncol=1, nrow=nrow(elec))
dta <- matrix(NA, nrow=1060, ncol=nrow(elec))
for(i in 1:nrow(elec))
{
x <- read.csv(as.character(elec[i,2]), header=F, skip=18, nrows=1060)
dta[,i] <- x[,2]
}
time<- x[,1]
#plot ladder and define peak threshold in HD
par(mfrow=c(3,4))
for(samples in 1:length(samples1))
#for(samples in 1:nrow(elec))
{
#which(dta[,lad] > 5)
#draw some segments and store max fluor in max.peaks matrix
lad<-samples
ti1<-17
ti2<-65
p.row<-length(dta[which(time==ti1):which(time==ti1+1),lad])*(ti2-ti1+1)
peaks <- matrix(NA, nrow=p.row, ncol=1)
max.peaks <- matrix(NA, nrow=length(seq(ti1,ti2,.05)), ncol=1)
for (i in 1:length(seq(ti1,ti2,.05)))
{
ii<-seq(ti1,ti2,.05)[i]
x0<- ii
x0<-round(x0,digits=3)
y0<- ii
x<- ii+.05
x<-round(x,digits=3)
y<- ii
#segments(x0,y0,x,y)
max.peaks[i]<- max(dta[which(time==x0):which(time==x),lad])
}
#human or mouse
cent<-max(max.peaks)*.6
#flatworm
#cent<-max(max.peaks)*.05
plot(time,dta[,lad],type="h",col="dodgerblue",xlab="seconds",ylab="Fluorescence Units",xlim=c(17,75),ylim=c(0,(max(max.peaks)*1.2)),main=elec[lad,1])
abline(h=cent,col="grey")
#define the middle of the segments
max.avg<-seq(ti1+.025,ti2+.025,.05)
#plot the maximum value of each segment
max.peaks
lines(max.avg,max.peaks,type="b",col="deepskyblue")
#lines(max.avg,min.peaks,type="b",col="blue")
#find time associated with peaks > 5 fu and fill seconds matrix with them
#12.5% of max fu incase chip is fu
p1<-which(max.avg==38)
p2<-which(max.avg==55)
sizes <- which(max.peaks>cent)
max.avg[sizes]
gt35<-which(max.avg[sizes]>35)
lt55<-which(max.avg[sizes]<55)
#as.numeric(duplicated(c(gt35,lt55)))
tfr<-which(duplicated(c(gt35,lt55))==TRUE)
c(gt35,lt55)[tfr]
max.avg[sizes][c(gt35,lt55)[tfr]]
max.peaks[sizes][c(gt35,lt55)[tfr]]
points(max.avg[sizes+1][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes-1][c(gt35,lt55)[tfr]],col="lawngreen",pch=19)
points(max.avg[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],col="deeppink3",pch=19)
#colnms2<-c("max.avg[sizes]","max.peaks[sizes]","max.avg[sizes+1]","max.peaks[sizes+1]","max.avg[sizes-1]")
max.avg[sizes]
max.peaks[sizes]
max.avg[sizes+1]
max.peaks[sizes+1]
max.avg[sizes-1]
max.peaks[sizes-1]
area<-cbind(max.peaks[sizes-1][c(gt35,lt55)[tfr]],max.peaks[sizes][c(gt35,lt55)[tfr]],max.peaks[sizes+1][c(gt35,lt55)[tfr]])
time2<-cbind(max.avg[sizes-1][c(gt35,lt55)[tfr]],max.avg[sizes][c(gt35,lt55)[tfr]],max.avg[sizes+1][c(gt35,lt55)[tfr]])
p1<-area[which(area[,1]<cent),1][1]
t1<-time2[which(area[,1]<cent),1][1]
p2<-area[which(area[,1]<cent),1][2]
t2<-time2[which(area[,1]<cent),1][2]
#area[which(area[,2]<cent),2]
#area[which(area[,2]<cent),2]
p3<-area[which(area[,3]<cent),3][1]
t3<-time2[which(area[,3]<cent),3][1]
p4<-area[which(area[,3]<cent),3][2]
t4<-time2[which(area[,3]<cent),3][2]
#abline(v=area[which(area[,3]<cent),3])
#abline(h=11)
#abline(v=42.75)
v1<-t1
v2<-t2
v3<-t3
v4<-t4
abline(v=v1)
abline(v=v2)
abline(v=v3)
abline(v=v4)
#time[area[which(area[,1]<cent),1][1]==dta[,lad]]
#time[area[which(area[,3]<cent),3][1]==dta[,lad]]
#time[area[which(area[,1]<cent),1][2]==dta[,lad]]
#time[area[which(area[,3]<cent),3][2]==dta[,lad]]
#abline(v=time[area[which(area[,1]<cent),1][1]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][1]==dta[,lad]])
#abline(v=time[area[which(area[,1]<cent),1][2]==dta[,lad]])
#abline(v=time[area[which(area[,3]<cent),3][2]==dta[,lad]])
#abline(h=area[which(area[,3]<cent),3][1])
#abline(h=area[which(area[,3]<cent),3][2]+3)
s18a<-round(t1,digits=1)
s18b<-round(t3,digits=1)
s28a<-round(t2,digits=1)
s28b<-round(t4,digits=1)
if(is.na(s18a)){s18a<-30}
if(is.na(s18b)){s18b<-35.05}
if(is.na(s28a)){s28a<-30}
if(is.na(s28b)){s28b<-30.05}
eta<-which(time==s18a)
etb<-which(time==s18b)
tea<-which(time==s28a)
teb<-which(time==s28b)
#dta[eta:etb,lad]
#dta[tea:teb,lad]
s18<-sum(dta[eta:etb,lad])
s28<-sum(dta[tea:teb,lad])
s28/s18
qc.mat[samples]<-  (-1*exp((s28/s18)*-1)+1)*10
text(60,cent*.75,c(format(qc.mat[samples],digits=2),"\n\nB-RIN"))
}
qc.mat
###end of ratio calc
###
###
#barplot(qc.mat,beside=T,ylim=c(0,11))
#set the working directory
#results <- read.delim("results-files.txt", sep=",",header=F,skip=0,as.is=T)
#targets <- read.delim(results[1,], sep=",",header=F,skip=14,as.is=T)
#targets[1:5,]
tbl<-which(targets[,1]=="Overall Results:")
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
#read in the targets file
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
mat1 <- matrix(NA, ncol=length(colnms), nrow=length(samples)*length(results))
colnames(mat1)<-colnms
rownames(mat1)<-rep(samples,length(results))
for(i in 1:length(results))
{
targets <- read.delim(results[i], sep=",",header=F,skip=14,as.is=T)
targets[1:5,]
which(targets[,1]=="Sample Name")
length(which(targets[,1]=="Sample Name"))
samples<-as.character(targets[which(targets[,1]=="Sample Name"),2])
tbl<-which(targets[,1]=="Overall Results:")
length(which(targets[,1]=="Overall Results:"))
tbl[1]
tbl[1]
colnms<-as.character(targets[(tbl[1]+1):(tbl[1]+4),1])
for (j in 1:length(samples))
{
mat1[((length(samples)*(i-1))+j),]<-as.matrix(targets[(tbl[j]+1):(tbl[j]+4),2])
}
}
mat1
cbind(mat1,qc.mat)
colnames(qc.mat)<-"B-RIN"
write.table(cbind(mat1,qc.mat),file=paste(dir,"\\RIN_VS_B-RIN.txt",sep=""), sep="\t",col.names=NA)

Thanks for posting and answering your own question! It would be great if you could also explain how to use this code to do what your question asks for.


text:
Building network for biomarker interaction analysis
I am looking for advice on how to best approach a problem I am faced with.
I have a dataset of numerous degradative biomarkers, clinical information and various other measures (from clinical trials). I would like to see how different biomarkers relate to each other, and effect their levels (high/low etc).
My goal is to gain an understanding of how these biomarker levels fluctuate in different disease profiles - and understand if there is a relationship between each of them (co/independent). The limitation is that I do not have access to genetic or environmental data so this will be a "random effect". There are around 1000 marker samples and 10's of biomarkers mixed between clinical and biochemical. There are also numerous time points to be incorporated.
I have thought about implementing something similar to a protein-protein interaction network, although this will of course not involve physical interactions, but more pathway interactions for further hypothesis driven research.
I wanted to firstly: validate if this project actually makes sense and secondly: ask for some advice around how to best implement this.
I have experience implementing neural/deep nets, directed and undirected networks (Bayes) in python, but not sure exactly how to go about this.
Thanks!

First of all it is important to know very well what kind of data you have and the limitations and advantages it has. It is also equally important to know your goal: which is your goal? Find the best biomarker? I was starting to answer, but then I thought that asking this won't damage my answer. Also which is the number of samples for each data type?

Yes agreed! My goal is to gain an understanding of how these biomarker levels fluctuate in different disease profiles - and understand if there is a relationship between each of them (co/independent). The limitation is that I do not have access to genetic or environmental data so this will be a "random effect". Its not necessarily "best" biomarker as different diseases will have different marker profiles :) there are around 1000 marker samples and 10s of biomarkers mixed between clinical and biochemical. There are also numerous time points to be incorporated.

https://www.reddit.com/r/bioinformatics/comments/6n0x2s/building_network_for_biomarker_interaction/   Someone alse answered here :)

It would be fine if you [edit] the post to include the information from the previous comment and from reddit in your post. It will improve the question and make it self contained

You can start by looking at the correlation (if data is non-parametric) of each variable with each variable of each data point you have. One thing that might change is the correlation between variables along the disease.
Next you can determine which relationships are constant along the disease progression/clinical trial. Identify also the ones that don't show a constant relationship. To do so you can use WGCNA (although with just 10 variables is too low to expect a scale free network you can use to correlate the variables and group them by correlation profile)
Then the question is which variables are more important, those who kept the same relationship or those which change?
Use the eigengenes of those variables (separately by those that have a constant relationship and those which don't) to estimate the dependency with the parametric variables. Or you simply can do an ANOVA of all the variables with each non-parametric variable.

Thanks for your help, I am perhaps slightly further along the road/ am looking at something slightly different. I have looked at correlations between all variables be they biomarkers or end points etc. What (for me) is important to look at are co-dependencies of biomarkers across the whole system - all be it only 10 markers. I am here in the process of literature research :) https://arxiv.org/pdf/1501.04731.pdf

If you read carefully my answer. correlations are the start of the process to find the relationship between markers and profiles

Ok i see what you mean. I will have to look at WGCNA to see if it is extendable outside of genetic data. I am happy to implement this myself. I primarily wanted to make sure what I was proposing wasnt stupid (have poele done thins before - yes ) and what others experience was. Thanks for the advice!

Yes, it is extendable outside genetic data

@JB1 Perhaps if you wait a bit before selecting my answer as the best one you will get more answers (I doubt it, but is in general recommended). BTW if you find other methods let me know I am also interested in this area and will be doing a similar project soon

thanks for the advice! Ill let you know how i get on, libraries i find/ build


text:
Is there a way to retrieve ENSEMBL IDs from a search query?
I would like to be able to batch download FASTA files from ENSEMBL. I normally would use the API to download them from the ENSEMBL gene IDs. However, the IDs aren't ones I know, only ones I would get from searching on the ENSEMBL website and collecting the IDs manually.
Is there a way I can use the API to query search terms and collect all results that are returned?
So for example, if I use the ENSEMBL website to search "HLA" I could collect the list of IDs (ENSG00000204252 etc. etc.) manually. HLA works okay using the API to query gene IDs because there's only a few terms (A, B, C plus non-classical for class I) but I was wondering if there was a way to directly access the search query one the homepage programmatically, as I may need to do it for some messier examples.

Can you give an example of the sort of search you're running? There's lots of things in the API that might do what you want but don't think there's a generic "search" call.

Yes, sorry I didn't think to do that. I want to search a term associated with a gene family, say "HLA" and retrieve ENSEMBL IDs (or the sequences associated with, but I already have some code to do that) of the results that come up from the search.

What kind of terms are you looking for: Diseases terms, protein family names (KRT7,, KRT19, KRT11), ...? Why do you need to do so (maybe we can work around this problem) ?

Please [edit] your question and give us a few example IDs and the output you would like to get from them.

Yes, this is possible but you really need to give us specific examples. You say you query for HLA, OK, but what do you want to getr from it? The FASTA sequence of what? The gene? Primary transcript? All transcripts? Peptide? All peptides? What species? All species? This will almost certainly be possible either directly from the API or, at worst, by processing the HTML but we can't help unless you tell us what you need.

Sorry I've been trying to keep it relatively agnostic. I just need the ENSEMBL gene IDs for the human alleles. The problem is not what happens downstream of that, I just wanted a way of getting IDs of genes and their alleles in one go rather than searching each indivdual gene and getting the gene ID.

Yes, but it's hard to understand what you need without an example. You say you want to download FASTA, but also that you want to retrieve IDs. If you give us an example, we can show you how to do it in a way that is applicable to anything but we need an example to understand exactly what you need.

I would download the complete geneset FASTA file from ensembl and then parse out the genes you want with a script.
Alternatively, you could do similarly with API. Get all the gene names, select the ones you want, and then pull down each sequence.


text:
What sequencing data metrics should I record?
I'm new to bioinformatics and Im starting a new microbial sequencing project and I'd like to record all of the qc data correctly. This is research and I'm not sure what I'll need later. Is there any consensus on useful qc metrics to record in a database for each sequencing run?
No rnaseq, running FastQC, assembling with a defined reference for a single species.

The metrics from FastQC are more of a snapshot. It depends on your sequencing platform used. Was it illumina, iontorrent, pacbio?

The sequencing platform is illumina. Primarily miseq. I just don't want to regret not recording information that might be useful later.

Well all the important info for illumina will be stored in basespace you would be looking at cluster density and the Q-score.

Is there any reason not to keep all of it?

No, just wanted to distill the essentials into a local database if possible.

What it might not be essential today might be essential in a couple of years. You can always delete it later, but you can't recover data after years of the experiment

If you are using a MiSeq is your data stored on basespace or locally?  If it is on basespace then the QC metrics should be store there for you.  If locally then you can go into the run data and pull the metrics files from the Illumina reports.  Really you are looking for information about clustering etc.  The sequencing quality and read metrics can be recalculated from the fastq/Bam files at any time.

@AaronBerlin we're storing locally. Thanks for the suggestions.

This also depends what you want to do with the data. The FASTQs already have all the quality scores stored, so it's easy to re-create all the QC outputs if necessary. MultiQC is good at making these reports: http://multiqc.info/
If you're planning to assemble them, the metrics matter a bit more. You'd want to record things like coverage depth, number of contigs, length of contigs etc. People usually report N50 for example

I would say that the following metrics are useful:

Sample source (Extraction place, person, date)
Details about pre-HTS preparation (e.g. amplicon PCR, rRNA removal)
Reads per sample (+ per barcode / lane, if available)
Sequencing kit (e.g. TruSeq v3), useful for adapter trimming
Sequence length (e.g. 125bp)
Paired end / single end
If paired end, target fragment length in sample prep (e.g. 650bp)
Reference sequence metadata

Scientific name
Location of reference genome (preferably a stable web address)
Version of reference genome
Read mapping rate (per sample / barcode / lane), as a percentage of total reads
Mapping program + version

If available, potential contaminating organisms

Name / description
Location of contaminant sequence (preferably a stable web address)
Estimated mapping rate (bearing in mind that there may be some overlap with the target)
Mapping program + version

More generally, ENA / EBI has an excellent reference in sample checklists, which indicate the minimum acceptable metadata for attaching to sequence archive data:
https://www.ebi.ac.uk/ena/submit/checklists


text:
What tools can I use for a bacterial core/pan genome pipeline?
I want to perform a genome comparison on a group of isolates. I want to look into two broad groups of taxa and compare the accessory genome in each group. I have been using prokka (v1.12) and roary (v3.8.2) to do this but it appears the accessory_binary_genes.fa file is actually an untrue representation.
Note: gene_presence_absence.Rtab does contain all the full presence/ absence for accessory gene sets. Despite this Im still unhappy with the nomenclature of the gene groups [issue for another day]

[github issue 335] Your best to ignore the accessory_binary_genes.fa file. It is just for creating a quick and dirty tree with FastTree. The file itself is filtered to remove very common and not common variation to speedup the tree generation, hence the difference in numbers.
The top 5% and bottom 5% are excluded. It is truncated at 4000 genes.

I've been looking into alternative pipelines and a new software BPGA looks promising. Does anyone have experience with this?
I essentially want a tool which will give me the core and accessory gene sets, without the noise from partial gene hits.

You could use OrthoMCL which is an older and less easier to use implementation of the same workflow. I believe that will give you all the gene names back rather than truncating anything.

You can still use roary and use the information contained in the gene_presence_absence.csv file, which is complete and contains the gene IDs for each orthologous group

LS-BSR should be able to give you what you're looking for. See the article.
After you run the primary analysis, there is a simple documented workflow for splitting the pangenome into core and accessory, based on a user defined threshold. I'm the developer, so can help if you run into any problems.

Hi Jason, Thanks for the answer. Isnt LS-BSR meant to be less sensitive than roary? Also could you provide a link/ source for the workflow you mention. Thanks.

Also how does LS-BSR address genes broken due to contig borders?

This will obviously cause problems. The truncated gene will appear to be divergent or missing depending on where it is broken. One way that we address this problem is to map the raw data back against the pan-genome and determine if the region is truly missing or not, although this is not part of the LS-BSR package

I think these are really important insights that you should integrate to your answer by [edit]ing it. Also, do not hesitate to use hyperlinks to your scripts, but try to give a complete answer, so next time somebody will have similar pan genome problem, he will understand what to do now.

Roary also takes into account paralogs, so sometimes two core genes are split into different groups based on their neighbour genes and they end up with different nomenclature (group_*...). As suggested by Andrew Page in the github issue I would consider the gene_presence_absence.Rtab (this contains all the orthologous genes) and remove rows corresponding to vectors only containing 1s (core genes). In this way you will have a matrix of 1 and 0 corresponding to presence/absence of a particular accessory gene in your isolates.

Hi Sergio, yes this correct but my main concern is with gene names. If a gene is broken/ split between contigs it could end up labeled as e.g. Gene_1 and Gene_2 this would then be a mis representation of the genetic content. its partially an issue with prokka too Id imagine.

You are right, contigs border can lead to a false representation of the gene content.

Apparently there could be a new software to address these issue coming soon https://github.com/tseemann/prokka/issues/244


text:
How to convert GFF3 to GTF2
I would like to convert a file in gff3 format to a gtf2.2 format.
The reason why I would like to do this is: I have a set of transcripts assembled by a bunch of different software (and using RNA-seq data from different sequencing technologies) and I would like to compare each of those sets with the reference transcriptome annotation for that species (D. melanogaster). I already asked for the community advice about how to proceed on that, but to run the suggested software (TACO and cuffmerge) I need to have a GTF file containing the experimental transcripts to compare to the reference transcripts (GTF).
Up to now, I unsuccessfully tried:

gffread utility in the Cufflinks package (gffread input.gff3 -T -o output.gtf): this results in an empty output.gtf file and an empty log file (used Cufflinks v.2.2.1) - I contacted the authors via their Google group but haven't heard of them yet
gff3_to_gtf utility in the GenomeTools (gt) package (gt gff3_to_gtf input.gff3 -o output.gtf): the output is not created and the log file is not informative - I contacted the authors via their mailing list but haven't heard of them yet
GFF3_to_GTF utility in the FML package (./gff3_to_gtf_converter.pl input.gff3 output.gtf): the output just contains a header (##gff-version 2.5) and the log is empty

The gff3 file was created as output of GMAP, and contains the transcripts as found by alignment to the reference (specifying option -f gff3_match_cdna).
[Edit: what I found out a poteriori, is that this format is not a standard gff3, thus the conversion is not trivial...]
Here is the beginning of the gff3 file I tried to convert:
$ head r9_gmap_match-cdna.gff
##gff-version   3
# Generated by GMAP version 2017-04-24 using call:  /home/aechchik/software/gmap-2017-04-24/bin/gmap.sse42 -d gmapidx -D /scratch/beegfs/monthly/aechchik/isoforms/ref/chromosomes/ -f gff3_match_cdna -n 0 -t 20 r9_2d.fasta
2L    gmapidx    cDNA_match    18442664    18443024    79    -    . ID=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1;Name=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d;Target=ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d 141 489;Gap=M13 I10 M8 D3 M8 D3 M32 D1 M4 D2 M5 D3 M34 I3 M6 D1 M7 D1 M5 D1 M30 D2 M16 I1 M8 D3 M10 D1 M5 D1 M6 I1 M8 I1 M8 D1 M33 D3 M33 I1 M28 D3 M25
###
3R    gmapidx    cDNA_match    15853880    15855465    96    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 80 1645;Gap=M46 D2 M12 D2 M157 D3 M4 D1 M50 I2 M3 I1 M12 I1 M68 I1 M66 I1 M53 D2 M47 D1 M16 D1 M35 D1 M155 D1 M28 D1 M166 D2 M47 D1 M8 D4 M69 D1 M28 D1 M5 D1 M12 D1 M202 I1 M115 D1 M61 I2 M7 D1 M7 I1 M36 D2 M41
3R    gmapidx    cDNA_match    15855529    15855742    97    +    . ID=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1;Name=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d;Target=dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d 1646 1856;Gap=M21 D1 M68 D1 M85 D1 M37
###
X    gmapidx    cDNA_match    14837810    14838142    93    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 74 406;Gap=M13 D1 M182 I1 M14 I2 M56 I1 M30 D2 M21 D1 M13
X    gmapidx    cDNA_match    14837470    14837753    92    -    . ID=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1;Name=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d;Target=960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d 407 688;Gap=I2 M5 I1 M64 I1 M44 D1 M9 D5 M31 D3 M23 I1 M19 I1 M25 I1 M26 D1 M19 I1 M9
###

I can't get cufflinks to work with this either. Strange. Could you explain what you need, exactly? I mean, GTF is basically GFF2. Would it be OK to just remove the attributes? What information do you need to keep? Basically, what will you use the GTF for? Do you need the gene_id and transcript_id fields and, if so, where would that information come from? It might be helpful if you came into our chat room. If you do, ping me (@terdon) to let me know.

I tried to convert this with rtracklayer in R and also failed (rtracklayer doesn't want to output proper GTF format) :( I wonder if a bit of python/perl would suffice.

The following bit of python code should work:
#!/usr/bin/env python
import sys

lastTranscript = [None, None, None, []]  # ID, chrom, strand, [(start, end, score), ...]

def getID(s):
    """Parse out the ID attribute"""
    s = s.split(";")
    for k in s:
        if k.startswith("ID="):
            return k[3:]
    return None

def dumpLastTranscript():
    """Print the last transcript"""
    bounds = sorted(lastTranscript[3])
    print("{}\tgmapidx\tgene\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    print("{}\tgmapidx\ttranscript\t{}\t{}\t.\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], bounds[0][0], bounds[-1][1], lastTranscript[2], lastTranscript[0], lastTranscript[0]))
    for start, end, score in bounds:
        print("{}\tgmapidx\texon\t{}\t{}\t{}\t{}\t.\tgene_id \"{}\"; transcript_id \"{}\";".format(lastTranscript[1], start, end, score, lastTranscript[2], lastTranscript[0], lastTranscript[0]))

def handleLine(cols):
    """Handle a single line, appending the exon bounds to the previous if relevant"""
    ID = getID(cols[8])
    assert(ID is not None)
    if lastTranscript[0] is not None and lastTranscript[0] != ID:
        dumpLastTranscript()
        lastTranscript[3] = []
    lastTranscript[0] = ID
    lastTranscript[1] = cols[0]
    lastTranscript[2] = cols[6]
    lastTranscript[3].append((int(cols[3]), int(cols[4]), cols[5]))

f = open(sys.argv[1])
for line in f:
    if line.startswith("#"):
        continue
    cols = line.strip().split("\t")
    handleLine(cols)

dumpLastTranscript()
f.close()

If you saved that as gff2gtf.py and made it executable, the usage would be ./gff2gtf.py foo.gff3 > foo.gtf. With the example that you provided in your post, the result is:
2L  gmapidx gene    18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx transcript  18442664    18443024    .   -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
2L  gmapidx exon    18442664    18443024    79  -   .   gene_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1"; transcript_id "ae6a7818-85b5-4739-8031-e58f4462ad41_Basecall_2D_2d.path1";
3R  gmapidx gene    15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx transcript  15853880    15855742    .   +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15853880    15855465    96  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
3R  gmapidx exon    15855529    15855742    97  +   .   gene_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1"; transcript_id "dd2444cf-34d6-4cd3-87ab-0ae1f3cb1f96_Basecall_2D_2d.path1";
X   gmapidx gene    14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx transcript  14837470    14838142    .   -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837470    14837753    92  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";
X   gmapidx exon    14837810    14838142    93  -   .   gene_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1"; transcript_id "960b50cd-945e-4c12-b9bc-367f965575bb_Basecall_2D_2d.path1";

That at least looks like correct gtf 2.2 format to me.

Close to be GTF2.2 but is not. Feature type from 3rd column should be filtered to match the GTF2.2 specification. Only 9 feature types are accepted. Gene is not as example. It became accepted in version 2.5. See here https://agat.readthedocs.io/en/latest/gxf.html#gtf2-2

@juke34 As there is no coherent GFF specification, you're always going to have to tweak something or other.

There are clear specifications (I would say excepted for GFF0 and GTF1) but it is easy to not follow them ^^. Have you look at the review I made? I tried to gather all the information in one place:  https://agat.readthedocs.io/en/latest/gxf.html

The GenomeTools package has lots of small applications that let you massage GFF3. It may have a converter, or it may let you tweak your existing GFF3 so that gffread works

Oh... I should read past the first item on the list... Try using gt gff3 to transform and output GFF3 files?

yes @Dan, I tried gt (Genome Tools package) - see point 2 of my question, and it is not working as it is.

Please edit your question and include an example of how the OP can use your suggested tool to do what they ask. As it stands, this is more of a comment than an answer. Please also see https://meta.stackexchange.com/q/225370/203101

@terdon Do you mean the answer?

@KamilSJaron d'oh! Yes, sorry, I was asking Dan to edit this answer, not the question.

I'm quite sure it will work with agat_convert_sp_gff2gtf.pl from AGAT.
See the review of the different solution I made here:
https://agat.readthedocs.io/en/latest/gff_to_gtf.html
AGAT it the only one you can decide wich version of GTF you require as output.


text:
Bioperl - how can i print first result of search sequence per iteration?
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;

$query = "LEGK";
$query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );

$gb_obj = Bio::DB::GenBank->new;

$stream_obj = $gb_obj->get_Stream_by_query($query_obj);

while ($seq_obj = $stream_obj->next_seq) {
    # do something with the sequence object
    print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq[,'\n';

How can I print first occurrence of protein sequence?

Please [edit] your question and clarify. Do you mean print just the first sequence? Or do you mean that if multiple identifiers have the same sequence, you only want to print the first of them? Or is it that you might have different identifiers with the same sequence and you only want the first of those?

In the simplest case, if you just want to stop after the first record was printed, you can just add exit (I also corrected the syntax errors you had and added use strict and use warnings; I suggest you get into the habit of using those two, they save you from a lot of grief in the long run):
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my $query = "LEGK";
my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
my $gb_obj = Bio::DB::GenBank->new;
my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
while (my $seq_obj = $stream_obj->next_seq) {
 print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
    exit;
}

If instead, you want to only process the first result for a list of multiple queries:
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my @queries = qw(LEGK TP53 ACT1);
query:foreach my $query (@queries) {
 my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
 my $gb_obj = Bio::DB::GenBank->new;
 my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
 while (my $seq_obj = $stream_obj->next_seq) {
  print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
        ## Move to the next query
        next query;
    }
}

Or, to print only the first sequence if multiple quries have identical sequences:
#!/usr/bin/env perl
use Bio::DB::GenBank;
use Bio::DB::Query::GenBank;
use strict;
use warnings;
my @queries = qw(LEGK TP53 ACT1);
foreach my $query (@queries) {
 my %seen;
 my $query_obj = Bio::DB::Query::GenBank->new(-db => 'protein',
                                          -query => $query );
 my $gb_obj = Bio::DB::GenBank->new;
 my $stream_obj = $gb_obj->get_Stream_by_query($query_obj);
 while (my $seq_obj = $stream_obj->next_seq) {
  if ($seen{$seq_obj->seq}) {
   ## Here, since we're not using an argument (next query), this will move
   ## to the next iteration of the current loop and not the `foreach`.
   next;
  }
  print ">$query",' ', $seq_obj->display_id, ' ', $seq_obj->desc,"\n", $seq_obj->seq, "\n";
    }
}


text:
RNAseq: Z score, Intensity, and Resources
I'm very new to bioinformatics in general, and I'm trying to understand some basic concepts.
I have RNAseq data, and bioinformatics people tell me that intensities cannot be compared across patients. So there are all of these pipelines to compare intensities to Z scores--are those as simple as just plugging data into a bioconductor package?
It would be great to get some overview description of why Z values are important/why you can't compare intensities across patients, and/or some pointers towards any resources that I can read to learn more.
Thanks much for your time!

Your wording suggests that you have microarray data rather than RNAseq, is this correct? Z values are rarely used in RNAseq, rather counts per gene/feature or TPMs are typically used.

I have some of both--even though I think I've seen Z-scores from TCGA illumna sets? Or perhaps I was mistaken.

It depends on what test or analysis you want to do, whether you need intensities (expression values) or z-scores.
If you want to do statistical analysis, such as finding differentially expressed genes between groups of patients (e.g., with limma), you don't want to use z-scores. But you use normalized intensities (for microarrays) or start with raw counts for RNAseq. Then follow the user guide (limma has a great user guide, also for starters).
For visualization in heatmaps or for other clustering (e.g., k-means, fuzzy) it is useful to use z-scores. Z-scores are a form of transformation (scaling), where every genes is sort of "reset" to the mean of all samples, using also the standard deviation. If you want to know exactly what a z-score is, a simple google search can tell you the details.
In R you can use the scale function for z-score transformation. Be aware that the function works on columns though. Which means you have to transpose your matrix first if genes are in the rows, and then transpose them back after transformation.

Rather than transposing the matrix I’d suggest using apply with scale.

In cancer genomics experiments often involve sequencing a large number of tumour samples, but few or no matched normals. This is partly due to financial constraints, but also due to ethical ones: it is easy to get ethical approval for tumour tissue, which is removed as a matter of course in a patient's care anyway, whereas it is hard to get approval for normal tissue, which would not normally be removed and would require an additional procedure.
However, this presents a problem. Most of the procedures we use in standard gene expression analysis involve comparing one group to another after suitable normalisations of the data. But here we have no suitable comparison group. To get around this problem in cancer RNA-seq analysis we do something called outlier detection. That is, we assume that for a given gene, the majority of the patients have "normal" expression, and look for the few patients that don't look like the rest. We do this with the Z-score.
To calculate the Z-score for an observation we subtract the mean of all observations and divide by the standard deviation. Thus the Z score of an observation is how many standard deviations an observation is from the mean of all observations - or how unusual it is.
Thus for Gene A in Patient 1, we calculate how many standard deviations Gene A is from the mean of Gene A across all patients. A very big (or small/negative) value tells us that the expression of Gene A is unusual in patient 1 compared to the other patients.
(NB: in some cancer studies they use the median and median absolute deviation rather than the mean and standard deviation to compute a "robust Z-score" as the data is unlikely to be normally distributed).

Interesting and well explained. Thanks. About the use of the median, I would add that this makes the results less sensitive to outliers.


text:
Chip-Seq data description
I have the following Chip-Seq data and could not found a description of it. Can you help me with it? More, I would like to find out the order of the nucleotides in binding regions; that is, what number they are.
@GWZHISEQ01:319:HA1NPADXX:1:1101:1182:2196 1:N:0:ACTGAT
GGAGGAAGTGTGTTGCTGCAGACAGTAGCAACCAGACCCACACTGCGCGTA
+
CCCFFFFFHHFHHIJJJJIJJIIIIHGJJJJJJJCGHIJJIJJFHIIJJII
@GWZHISEQ01:319:HA1NPADXX:1:1101:1272:2109 1:N:0:ACTGAT
NCCACGCTAGGCTCAGCTTGTCGGCCTGGCTAAGCAGTTGCGAAAGTGCGC
+
#1=DDDFFHHHHHJJIJJJJJJIGEHGGJIGHJFHGJIJJJJJJJ@CHCBG
@GWZHISEQ01:319:HA1NPADXX:1:1101:1418:2141 1:N:0:ACTGAT
GCACGCACTACCCAGAGATCATCCAAAGCCTGAAGCCACAGGGCGCACTCG
+
CCCFFFFFGHGHHJJJIIJJJJIIHGGGIJIIJJIIIJJJIJJJIJIJIJH
@GWZHISEQ01:319:HA1NPADXX:1:1101:1302:2160 1:N:0:ACTGAT
ACACTCTTTCCCTACACGACGCTCTTCCGAGATCGGAAGAGCACACGTCTG
+
CCCFFFFFHDHHHJJJJIJJJJIIIIIJIHJIJJJIIJHIIJJJJJJJJHH

Cross-posted on biostars

Is that a problem?

It's frowned upon at the very least.

I would argue that cross-posting is not a problem provided you reference the cross-posts on each site.

That's a fastq file, you will want to align it to the genome, call peaks, and then use something like MEME to determine binding motifs.

how do I align it to the genome?

@frosh You need more help than a Q&A site can offer. I recommend that you attend a workshop on NGS data analysis, or find a collaborator who can do the analysis for you.

You're right. I will delete the question.

You don't need to delete it.

The following mini-courses may be useful: http://rafalab.github.io/pages/harvardx.html


text:
Subset FASTA file by species name
I have a problem: I've managed to download a massive fasta file of 1500 sequences, but now I want to split them into separate fasta files based on the genus.
EDIT
The fasta file looks like this:
terminase_large.fasta
>YP_009300697.1 terminase large subunit [Arthrobacter phage Mudcat]
MGLSNTATPLYYGQF...

>YP_009208724.1 hypothetical protein ADP65_00072 [Achromobacter phage phiAxp-3]
MSNVLLKQELDEWL...

>YP_009148449.1 large terminase subunit [Delftia phage RG-2014]
MSEPRKLVKKTLD...
...

I would like to end up with something like this:
mycobacterium_terminase.fasta
>Mycobacterium phage JAMaL
MVRKKPPPELE...

>Mycobacterium phage Bruin
MEVCGYTLDDI...

>Mycobacterium phage Zaka
MSLDNHLPELA...

salmonella_terminase.fasta
>Salmonella virus SETP7
MNVDITATEPQ...

>Salmonella virus SE2
MEGGSDSLIAM...
and so on...

So that the label on the alignment would have the name of the phage and not the protein's ID.
I've managed to extract the genus names of my organism with this:
outfile = open('species.txt', 'w')
with open('terminase_large.fasta') as fd:
    for line in fd:
        if line.startswith('>'):
            if '[' in line:
                name=line.split('[')[-1]
                name=name.split(' ', 1)[0]
                outfile.write(name[:] + "\n")
outfile.close()

And got to extract only the unique names with this:
lines_seen = set()
outfile = open('species2.txt', "w")
for line in open("species.txt", "r"):
    if line not in lines_seen:  # not a duplicate
        outfile.write(line)
        lines_seen.add(line)
outfile.close()

(Can I merge those two scripts together?)
Now, my genus names look like this:
Arthrobacter
Achromobacter
Delftia
....

I tried automating my script to get the Entrez data, but it gives me the 'Supplied id parameter is empty' message
My code looks like this:
from Bio import Entrez
Entrez.email = "example@example.org"

for line in open("species2.txt", "r"):
    searchterm = "(terminase large subunit AND viruses[Organism]) AND" +line+ "AND refseq[Filter]"
    searchResultHandle = Entrez.esearch(db="protein", term=searchterm, 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('terminase_large_'+str(line[:-1])+'.fasta', 'w')
    out_handle.write(record.rstrip('\n'))

Can someone help me with it?

Please [edit] your question and i) show us a short example of your input file, including different genera and ii) show us the output you would expect from that file. That way, we can test our solutions and be sure to give you what you need.

Also, you seem to have bacterial genera, why are you searching for viruses? Of course that won't find any results.

My keywords use the bacteria genera and the virus as organism, since I want specific phages that infect bacteria, eg Enterobacter phage or Vibrio phage

Ah yes, sorry, I thought you were specifying the bacterial genus as [ORGANISM] as well. So, is it enough to parse the existing file or do you need your script to be modified to create the multiple files next time it runs? Modifying the existing one and splitting it into subfiles by genus is very easy. Would that be enough?

@terdon sure, I just didn't know which was easier: splitting the existing file, which has the genera mixed up or creating a new fasta file by searching Entrez

Your output and input don't match, I can't see "Mycobacterium" in the input file. Is this intentional?

@gringer yeah, it's confusing. The OP wants to i) change the fasta ID lines and ii) split into multiple files, collecting all sequences from phages infecting a given bacterial genus together. They only showed the output for the first part and not the second.

Sorry for the bad files, I will try to edit it once more

Thanks @tahunami, the output seems to fit better with your question now. Are you able to change the input so that it only includes information in the output file? I still can't see any of the output proteins or names in the input example. Feel free to put dots in to demonstrate "skipping over thousands of lines", but I think it's important that at least everything in the output example should be in the input example.

I'm not sure this solves your problem, but I would use line.strip() in searchterm, not just line.

Splitting into multiple files and changing the IDs can be easily done:
perl -pe 'if(/>/){/\[(.*?)\]\s*$/; $_="> $1\n"}' file.fa |
    awk '(/^>/){name=$2} {print >> name".fa"}'

That assumes all your FASTA headers have [foo bar baz] as the last element of a line. It will create a file called foo.fa (the bacterium's name) with all sequences saved there.
Explanation

perl -pe : run the script given by -e on each line of the input file, and print the resulting line.
if(/>/) : if this line starts with a >. \
/\[(.*?)\]\s*$/ : match an opening bracket (\[), then capture (that's what the parentheses do, they capture a pattern so we can refer to it as $1) everything until the first ] (.*?\])
$_="> $1\n" : the $_ special variable in Perl is (in this case) the current line. So, $_=foo means "make the current line read foo. Since the -p prints each input line, changing the value of $_ means the changed value will be printed. So here, we are printing >, whatever was in the square brackets ($1) and a newline character.

The output of the perl command alone on your example input file is:
$ perl -pe 'if(/>/){/\[(.*)\]\s*$/; $_="> $1\n"}' file.fa
> Arthrobacter phage Mudcat
MGLSNTATPLYYGQF...

> Achromobacter phage phiAxp-3
MSNVLLKQELDEWL...

> Delftia phage RG-2014
MSEPRKLVKKTLD...

So, we now pass it through awk which does:

(/^>/){a=$2} : if this line starts with an >, save the second field (the bacterial species) as the variable name.
{print >> name".fa"} : print each line into a file whose name is the current value of the variable name with a .fa. extension.

If you prefer python scripts to the one-liner approach, you can do the same thing with:
#!/usr/bin/env python
import re
outFile = None
with open("file.fa", "r") as inFile, open(species + ".fa", 'a') as outFile:
    for line in inFile:
        line = line.rstrip()
        if line.startswith('>'):
            regex = re.compile('.*\[((.+?)\s+.*?)\].*')
            matches = regex.search(line)
            species = matches[2]
            outFile.write('>%s\n' % matches[1])
        else:
            outFile.write("%s\n" % line)

As for your script, you've got the right idea, but have a small bug. You forgot to remove the \n from your input file, so it looks for Arthrobacter\n instead of Arthrobacter. The golden rule of debugging is "print all the things". If you add print("Searchterm: ",searchterm)  you will see:
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDArthrobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDAchromobacter
AND refseq[Filter]
Searchterm:  (terminase large subunit AND viruses[Organism]) ANDDelftia
AND refseq[Filter]

So, you need to remove the newline characters and add a space like so (I also made it a bit more "pythonic" and conforming to the Python syntax guidelines):
#!/usr/bin/env python
from Bio import Entrez
Entrez.email = "example@example.org"

with open("species2.txt", "r") as in_handle:
    for line in in_handle:
        line = line.rstrip()
        searchterm = ("(terminase large subunit AND viruses[Organism]) " +
                      "AND %s AND refseq[Filter]" % line)
        print("Searchterm: ", searchterm)
        searchResultHandle = Entrez.esearch(db="protein",
                                            term=searchterm, retmax=1000)
        searchResult = Entrez.read(searchResultHandle)
        ids = searchResult["IdList"]

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

        with open('terminase_large_' + str(line[:-1]) + '.fasta', 'w') as out_handle:
            out_handle.write(record.rstrip('\n'))

Is it possible to do the perl script in python? I don't have perl on my pc installed.
About the second part: it works! Sort of. When I debug my script I can see that python gets the desired fasta file (a url link is created), and it creates a file with the correct name, but somehow doesn't input the data into it

@tahunami what do you mean, unless you're using Windows, it should be installed by default. Are you sure it's not installed?

I've saved in itno a .pl file, and run it on perl, but got "Execution of split.pl aborted due to compilation errors"

@tahunami no, that is a one-liner. Don't save as a file, just paste it directly into an open terminal. That's just a *nix command. In any case, I also added a python version for you.

A (bio)awk-based solution
Generate the list of genus
bioawk -c fastx '{print $comment}' terminase_large.fasta \
    | sed -r 's/.*\[(\w+).*/\1/' \
    | sort -u > genus_names.txt

Here is how it works:

bioawk is like awk but with extra parsing capabilities. -c fastx parses the fields in fasta or fastq records. {print $comment} will print the $comment field which contain the genus names you want to extract.
The sed command captures (with parentheses) the word (\w+) that is just after the opening square bracket (which should be the genus name, if there is only one opening square bracket on the line), and substitutes the whole line (that is, everything before the bracket (.*), the bracket (\[), the word, and everything after the word (.* again)) with just the captured word (\1).
All this is sorted, and the -u option of sort ensures there is only one occurrence of each genus name in the output, which we put in the genus_names.txt file.

Extract the fasta records for each genus
for genus in $(cat genus_names.txt)
do
    bioawk -c fastx '{print}' terminase_large.fasta \
        | grep ${genus} \
        | awk -F "\t" '{print ">"$4"\n"$2}' \
        | sed -r 's/^>[^\[]+\[(.*)\]/>\1/' > ${genus}.fasta
done

Here is how it works:

We loop on the content of the previously established list of genus names.
For each value of ${genus}, we parse again the fasta file with bioawk. This time, for a given fasta record, we print all the elements that have been parsed. They appear on one line, in the following order (according to bioawk -c help): 1:name 2:seq 3:qual 4:comment (we actually only need seq and comment for the following steps, but there are less risks of making errors at the later awk step if we keep the four fields).
We use grep to select the resulting lines that contain the genus name we're currently dealing with.
The lines that have been successfully selected are then parsed by awk and reformatted into fasta format using only the comment part in the header. The -F "\t" is to indicate that the field delimiters are tabulations (which is how bioawk has written the records when we did {print}).
We use sed to keep only the part inside the square brackets (again using a pattern capture between parentheses, that we use in the substitution part with \1), and the output is written in a file named using the current value of ${genus}.

More details on the last sed command:
We match the whole header lines. They start with ">" (^>), then there are some non-"[" characters ([^\[]+), then a "[" (\[), then something we capture ((.*)), then a "]" (\]). We replace them with the captured part. The sequence lines are unaffected because they don't match the above regular expression.

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