text:
Curl error when installing Bioconductor packages with Microsoft R Open 3.5.0
I'm trying to install Bioconductor packages with Microsoft R Open 3.5.0 seems to fail not sure what is the issue.
This is what I'm trying to do
source("https://bioconductor.org/biocLite.R")
biocLite("BiocInstaller")

But I keep getting an error about "curl" had non-zero exit.
How can I install the packages of Bioconductor?

Please [edit] your question and clarify. Your first paragraph seems to be about installing R but the last one suggests you have already installed it. Please don't ask multiple questions in a single post, [edit] and ask one, specific and clear question. Thanks!

updated my question ..iim getting this "curl' non-zero exit which wanted to post the error but closed my terminal so lost the exact error message

Then please run the command again and then add the error message to the your question. We can't help you solve an error if you don't show it.

I can't comment yet but from what you said it seems that you need to install additional libraries like the one you mentioned: "curl".
from the R terminal you should type:
install.packages("curl")

then R will install the package. Not sure if there are some other dependencies. If yes R will return you some error message in which there are some package names you need to install.
EDIT:
of course, after installing all the required packages you should be able to install biocLite too.

cd /opt/microsoft/ropen/3.5.1/lib64/R/library/  inside the path i can see the curl so why still i get the error i'm not sure

try loading the curl package from the R terminal and see if it actually works (library(curl)). Maybe some libraries are missing. I use linux and curl requires additional libraries that are not installed by R but are system's libraries.

yeah i did check the system level libraries as i can see using this "locate libcurl "so still not sure why the non-zero exit

finally after some time it did work.

sorry for not replying! if you solved the issue just post how you did it. This will help other people having the same problems!

earlier it was not showing any error no idea what was going on but later it showed up the dependencies  error which i installed it sort of fixed the problem ..but the way i used to update using CRAN is now no longer working with Microsoft R open here i can;t define the repository how to fix that any suggestion ?

@krushnachChandra, sorry for this very late reply, I must have lost your comment...have you solved your issue? otherwise I believe you should open another question with more details.

well thank you for your late response , i had few other issues as well .So i reverted to the base R


text:
Splitting fasta file into smaller files based on header pattern
I have to split this fasta files into smaller files and write them into individual files my files
>lcl|CP000522.1_prot_ABO13860.1_1 [locus_tag=A1S_3471] [protein=hypothetical protein] [protein_id=ABO13860.1] [location=1..957] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13850.1_2 [locus_tag=A1S_3461] [protein=DNA replication protein] [protein_id=ABO13850.1] [location=950..1504] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13851.1_3 [locus_tag=A1S_3462] [protein=hypothetical protein] [protein_id=ABO13851.1] [location=complement(2523..3437)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13852.1_4 [locus_tag=A1S_3463] [protein=YPPCP.09C-like protein] [protein_id=ABO13852.1] [location=3538..4788] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13853.1_5 [locus_tag=A1S_3464] [protein=Cro-like protein] [protein_id=ABO13853.1] [location=5039..5629] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13854.1_6 [locus_tag=A1S_3465] [protein=hypothetical protein] [protein_id=ABO13854.1] [location=complement(6340..6906)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13855.1_7 [locus_tag=A1S_3466] [protein=Resolvase] [protein_id=ABO13855.1] [location=complement(7074..7685)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13856.1_8 [locus_tag=A1S_3467] [protein=hypothetical protein] [protein_id=ABO13856.1] [location=complement(8602..9732)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13857.1_9 [locus_tag=A1S_3468] [protein=putative lipoprotein] [protein_id=ABO13857.1] [location=complement(10072..10374)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13858.1_10 [locus_tag=A1S_3469] [protein=Diaminopimelate decarboxylase] [protein_id=ABO13858.1] [location=complement(10367..10723)] [gbkey=CDS]
>lcl|CP000522.1_prot_ABO13859.1_11 [locus_tag=A1S_3470] [protein=regulatory protein LysR] [protein_id=ABO13859.1] [location=complement(12076..12444)] [gbkey=CDS]

The other pattern is
>lcl|CP000523.1_prot_ABO13861.1_1 [locus_tag=A1S_3472] [protein=DNA replication protein] [protein_id=ABO13861.1] [location=1..951] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13862.1_2 [locus_tag=A1S_3473] [protein=hypothetical protein] [protein_id=ABO13862.1] [location=3048..4262] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13863.1_3 [locus_tag=A1S_3474] [protein=hypothetical protein] [protein_id=ABO13863.1] [location=4357..5133] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13864.1_4 [locus_tag=A1S_3475] [protein=hypothetical protein] [protein_id=ABO13864.1] [location=6197..8608] [gbkey=CDS]
>lcl|CP000523.1_prot_ABO13865.1_5 [locus_tag=A1S_3476] [protein=secretory lipase] [protein_id=ABO13865.1] [location=8705..9403] [gbkey=CDS]

So now my idea is how do i parse and write them into individual files such that CP000522 output written to one file and CP000523 written to another file so forth and so on.
So far what i understand is i have to match the pattern after >lcl
so there are other patterns like "LN997847" in the file
Not sure how to proceed tried it in R but failed
it can be done with awk and sed which i tried but i can;t define something that parse all header like takes into account CP as well as LN .
Any suggestion or help would be highly appreciated .
here is my file

This is very easy with a bit of scripting one script will do it all. We do have R experts here (best solution)....Otherwise I'll post the code. How do you want each file named (important no matter which approach you take)?

Here's a simple awk approach:
awk '{if(/^>/){split($1,a,"[|.]")}print >> a[2]".fa"}' Protein_FASTA.txt

Or, more concisely, just:
awk '/^>/{split($1,a,"[|.]")}{print >> a[2]".fa"}' Protein_FASTA.txt

When run on the file linked to in your question, that results in the following files:
$ ls
AP014650.fa  CP003848.fa  CP007713.fa  CP012005.fa  CP015122.fa  CP017645.fa  CP018422.fa  CP020594.fa  CP023021.fa  CP024577.fa  CP026712.fa  CP027245.fa  CP030108.fa
CP000522.fa  CP003850.fa  CP007714.fa  CP012007.fa  CP015365.fa  CP017647.fa  CP018678.fa  CP020596.fa  CP023023.fa  CP024578.fa  CP026748.fa  CP027529.fa  CP030109.fa
CP000523.fa  CP003887.fa  CP008707.fa  CP012008.fa  CP015366.fa  CP017649.fa  CP018679.fa  CP021322.fa  CP023024.fa  CP025267.fa  CP026749.fa  CP027531.fa  CU459137.fa
CP000864.fa  CP003888.fa  CP008708.fa  CP012953.fa  CP015484.fa  CP017651.fa  CP019218.fa  CP021327.fa  CP023025.fa  CP026126.fa  CP027121.fa  CP027532.fa  CU459138.fa
CP000865.fa  CP003907.fa  CP008709.fa  CP012954.fa  CP015485.fa  CP017653.fa  CP020573.fa  CP021348.fa  CP023027.fa  CP026127.fa  CP027122.fa  CP027608.fa  CU459139.fa
CP001183.fa  CP003908.fa  CP008850.fa  CP012955.fa  CP015486.fa  CP017655.fa  CP020575.fa  CP021783.fa  CP023028.fa  CP026128.fa  CP027124.fa  CP027609.fa  CU459140.fa
CP001922.fa  CP003968.fa  CP008851.fa  CP012956.fa  CP016296.fa  CP017657.fa  CP020576.fa  CP021784.fa  CP023030.fa  CP026129.fa  CP027179.fa  CP027610.fa  JN377410.fa
CP001923.fa  CP004359.fa  CP010398.fa  CP013925.fa  CP016297.fa  CP018144.fa  CP020577.fa  CP021785.fa  CP023032.fa  CP026339.fa  CP027180.fa  CP029570.fa  LN865144.fa
CP001938.fa  CP006769.fa  CP010399.fa  CP014216.fa  CP016299.fa  CP018255.fa  CP020580.fa  CP021786.fa  CP023033.fa  CP026340.fa  CP027181.fa  CP029571.fa  LN997847.fa
CP002523.fa  CP007578.fa  CP010400.fa  CP014217.fa  CP016301.fa  CP018257.fa  CP020585.fa  CP021787.fa  CP023035.fa  CP026705.fa  CP027182.fa  CP029572.fa  LT594096.fa
CP002524.fa  CP007579.fa  CP010780.fa  CP014292.fa  CP016302.fa  CP018333.fa  CP020589.fa  CP022284.fa  CP024125.fa  CP026706.fa  CP027243.fa  CP029573.fa  Protein_FASTA.txt
CP003501.fa  CP007580.fa  CP010782.fa  CP014293.fa  CP017643.fa  CP018334.fa  CP020593.fa  CP022285.fa  CP024419.fa  CP026708.fa  CP027244.fa  CP030107.fa

Explanation

if(/^>/){split($1,a,"[|.]") : if this line starts with a >, split the first field on any occurrence of either | or . and save the results in the array a. Since your header lines all start with >lcl|, then the string you are looking for and a ., this means that the second value in the a array will be your target string.
print >> a[2]".fa" : print (append, >>) the current line to a file called "whatever the name of this sequence is" (a[2]) and .fa. This is run for every line in your input file. Note that if you run the same command again, you will need to first delete the files created the first time. If you don't, because I am using the >>, you will just append to the existing files.

I thought about perl -e  one-liner ... but I don't think a multiple file output is doable. The awk one-liner is cool.

Personally, I would have used a "_" along with "." and "|", but who am I. The small advantage of a script is that if its rerun it will clobber the old run output, without fear of appending.

@MichaelG. you can do it easily enough in perl, it's just slightly more cumbersome since you need to explicitly open and close each file: perl -ne 'if(/>.+?\|(.*?)\./){$name=$1; }open(my $fh, ">>","$name.fa"); print $fh "$_"; close($fh)' Protein_FASTA.txt. I didn't include the _ because the OP seemed to want the text before the version (the .N) to be the file name. I'd have used the _ myself too.

Thanks! I had no idea a file handle could be opened from "perl -e" ... I prefer awk code for this one.

wow so many solutions...let me run this..

thanks for the detailed explanation each line of the code..really appreciate it

You're welcome :). @MichaelG. note that I had a mistake in the first version of that perl one-liner. You need >>, not > (I've edited the comment) and it still has the same issue with existing files. The only way around that is to store everything in memory and only write at the end: perl -ne 'if(/>.+?\|(.*?)\./){$name=$1;}push @{$k{$name}},$_; END{for $name (keys(%k)){open(my $fh, ">","$name.fa"); print $fh @{$k{$name}}; close($fh)}}' Protein_FASTA.txt

I suggest replacing the if (…){…} block with the more awk-ish /^>/ {split($1, a, "[|.]")}, which simplified the actual print look to just {print >> a[2] ".fa"}.

@KonradRudolph that's how I would do it, yes. I used this version since I thought it would be easier for the OP to understand since it's more explicit and only 4 characters longer. But OK, I added the shorter version too.

@terdon I commented because I find the shorter version more readable (and I would argue that it is objectively so). Admittedly it does require knowing a bit about awk but that shouldn’t be an impediment.

@krushnachChandra ... its a culture. Make a code, condense a code and be more perlish, awkish or pythonic (as above).

I agree with @terdon flow control within one-line is minimised outside awk. This  awkish-ness so takes time to grasp.

source: Stack Exchange

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

sort: votes


text:
removing mitochondrial read and unassembled "random" from multiple bam files
This is the one i used for a single bam file to filter its mitochondrial as well as unassembled read
samtools view -h SRR2920466.bam | awk '{if($3!= "chrM" && $3 != "chrUn"){print $0}}' | samtools view -Shb - > this.filter.bam

Is there any way i can do it on multiple files with parallel as i do to index and sort bam files using samtools .
So i tried this ,which is of course not working as its not able to pipe properly because it shows me error regarding additional fonts
ls *.bam | parallel "samtools view -h {} | awk '{if($3 != "chrM" && $3 != "chrUn"){print $0}}'| samtools view -Shb {.}.bam"

Any help with my above codes? As i get it there is wrong with my pipe .
After updating the samtools to version 1.9 i ran the small script
i get this error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                    ^ unterminated regexp
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:                                                   ^ unexpected newline or end of string

I have 30 bam files i tried with 2 to see if its working or not , after the suggestion of @WouterDeCoster i tried the above one after updating samtools . the basic samtool with egrep works but for multiple file i couldn;t do it.

Your awk command is inside double quotes ("") which means that all of the awk field variables (e.g. $3) are being interpreted by the shell. Since $3 is not set in the shell, your awk is trying to compare an empty string and that's why you get the error:
awk: cmd. line:1: {if( != chrM &&  != chrUn){print ./index_sort.sh}}
awk: cmd. line:1:      ^ syntax error

Note how the error message shows if( != chrM and not if($3 != chrM. Since you need to use double quotes here because you need both an outer (for parallel) and an inner (for awk) quote, and single quotes can't be nested, you will have to escape the awk variables:
ls *.bam | parallel "samtools view -h {} |
    awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb {.}.bam"

looks like i need to clear my fundamental  awk concepts let me run and will update you. now i ran your code i get output like this "�&Ct�wh���x�Y#��J�f���{/.i��7��%M�������C�̛�<kbb�������B��f����<#?���
:P�y���
u��O��r.�_2ؑ�xu�zl�,xϭ"2��Uo�ȩ�@ayi���
/0���/���۶�
�>�Xh� ��Ί�5h�f~���euSa1����;5��{���R�3L����p%��[<�1)M�y!�w��#���_�`�V5rj��v�}�
S-���y)�Nμ4��Jӕ�D���n�1��ޥ�g�����y���䒋�N��)���<:���mʲ�溎[�ʘ���K��[punit@localhost chrm]$ 62;c62;c62;c62;c^C
"

@krushnachChandra well yes, of course. You are printing a bam file which is compressed binary data. What were you expecting to see? If you want sam output, remove the b from the options you give to samtools view. The b makes it print bam output.

so i idea is after filtering the mitochondrial read and unassembled one i want  bam output written to file with bam format which will be using for downstream analysis.I don't want to print it on the terminal . tried it with -o to write the output seems like i couldn;t make it work..

@krushnachChandra -o should work. Also, you could redirect in the parallel command: ls *.bam | parallel "samtools view -h {} | awk '{if(\$3 != "chrM" && \$3 != "chrUn"){print \$0}}' | samtools view -Shb  > {}.bam". But as I said under your question, the new error is completely unrelated. You probably have a malformed bam file. But please ask a new question about that. Even better, run for bam in *bam; do samtools view -h "$bam" | head -n1 || echo "$bam file is corrupted"; done That should help you find the bad file.

will run it and see

its working ,the output is written to bam file

i have accepted the answer sorry i missed it to mark since i saw this thread active again .

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


text:
Is it possible to perform MinION sequencing offline?
I vaguely remember, that the original plan of Oxford Nanopore was to provide cheap sequencers (MinION), but charge for base-calling. For that reason the base-calling was performed in the cloud, and the plan was to make it commercial once the technology is established.
Of course, this limits the potential uses of MinION in the field, since huge areas do not have decent internet connection. Also, not all the data can be legally transferred to a third-party company in the clinical studies.
For example for the Ebola paper, they had to create a special version of their software:

An offline-capable version of MinKNOW, with internet ‘ping’ disabled
  and online updates disabled was made available to us by Oxford
  Nanopore Technologies specifically for the project

There are couple of third-party base-callers available today. I am aware of Nanocall and DeepNano, but since they are not official, it can be hard for them to keep-up with the latest versions of sequencers and cells.

Is it possible as of today to sequence offline without a special arrangement (like the Ebola one).
If not, what's the policy of Oxford Nanopore toward third-party base-callers? Are they going to help them, or to sue them eventually?

Short answer: yes, but you need to get permission (and modified software) from ONT before doing that.
... but that doesn't tell the whole story. This question has the potential to be very confusing, and that's through no fault of the questioner. The issue is that for the MinION, sequencing (or more specifically, generating the raw data in the form of an electrical signal trace) is distinct and separable from base calling. Many other sequencers also have distinct raw data and base-calling phases, but they're not democratised to the degree they are on the MinION.
The "sequencing" part of MinION sequencing is carried out by ONT software, namely MinKNOW. As explained to me during PoreCampAU 2017, when the MinION is initially plugged into a computer it is missing the firmware necessary to carry out the sequencing. The most recent version of this firmware is usually downloaded at the start of a sequencing run by sending a request to ONT servers. In the usual case, you can't do sequencing without being able to access those servers, and you can't do sequencing without ONT knowing about it. However, ONT acknowledge that there are people out there who won't have Internet access when sequencing (e.g. sequencing Ebola in Africa, or metagenomic sequencing in the middle of the ocean), and an email to <support@nanoporetech.com> with reasons is likely to result in a quick software fix to the local sequencing problem.
Once the raw signals are acquired, the "base-calling" part of MinION sequencing can be done anywhere. The ONT-maintained basecaller is Albacore, and this will get the first model updates whenever the sequencing technology is changed (which happens a lot). Albacore is a local basecaller which can be obtained from ONT by browsing through their community pages (available to anyone who has a MinION); ONT switched to only allowing people to do basecalling locally in about April 2017, after establishing that using AWS servers was just too expensive. Albacore is open source and free-as-in-beer, but has a restrictive licensing agreement which limits the distribution (and modification) of the program. However, Albacore is not the only available basecaller. ONT provide a FOSS basecaller called nanonet. It's a little bit behind Albacore on technology, but ONT have said that all useful Albacore changes will eventually propagate through to nanonet. There is another non-ONT basecaller that I'm aware of which uses a neural network for basecalling: deepnano. Other basecallers exist, each varying distances away technology-wise, and I expect that more will appear in the future as the technology stabilises and more change-resistant computer scientists get in on the act.
Edit: ONT has just pulled back the curtain on their basecalling software; all the repositories that I've looked at so far (except for the Cliveome) have been released under the Mozilla Public License (free and open source, with some conditions and limitations). Included in that software repository is Scrappie, which is their testing / bleeding-edge version of Albacore.

In short, yes offline basecalling is supported.
Oxford nanopore has made available an official offline basecaller, called albacore. If you have an account to the nanpore community website you can download it from here.
Concerning running it in the field, last time I checked live basecalling was still requiring a lot of computing power, so you might want to run the MinION in the field, but basecall once back a your lab / computing facility.
I'm not aware of any plans to sue the third-party basecallers, I think Oxford Nanopore is pretty positive about them :-)

Running the MinION in the field and base-calling back in the lab defeats the whole (well, part of) the purpose: namely, to get results instantaneously.

True in some cases, but sometimes you want to sequence samples in remote areas / places where it'd be impractical / impossible to go back to a lab with good quality DNA samples (i.e difficult access to dry ice for transport). The MinION makes it possible to sequence in the field, then analyse back in the lab.
EDIT: Also, for small sequencing run you can probably basecall with a good laptop (8+ cpus)

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


text:
Heteroplasmy Calculation
For calculating heteroplasmy in mtDNA I was using MitoSeek, the tool provides position wise heteroplasmy percentage. So, will it be wise to add up all sites to get a position-wise percentage in a single sample, or calculate it by adding up all counts (variants) from an aligned file and divide it by total read count?

I haven't touched this subject, hence this is a comment, but for me it is better the second way, However it is much different between method 1 (I think you might want to use the mean instead of add up the percentage) and 2? What do other papers in this field?

I wouldn't think either of the options are appropriate. Variants at different sites could be arising from the same mtDNA copy, so I presume one would need to do some sort of LD modelling or something like that to account for this (or accept a relatively crude estimate).

Heteroplasmy at any site (as long as it's a reliable variant) suggests that there is mitochondrial heteroplasmy, it doesn't matter where it is. In that case, the "mitochondrial heteroplasmy" statistic would simply be the maximum value at any [reliable] location.
Due to the low amount of variation in the mitochondrial genome, and because it's a haploid genome, I don't think that a calculation like "average heteroplasmy" makes much sense for mitochondria.


text:
Understanding the significance of BLAT score calculations for read fragment alignment
I'm attempting to reconstruct the read fragments that form a hit to a target, excluding/skipping gaps. BLAST's outfmt 6, and alternatives such as Diamond's M8 format only give the number of gap opens on the alignment. One can't just use the start and end positions of the alignment as indices to slice the read DNA sequence, as this will include the nucleotides that were gapped by the alignment.
PSL appeared to be the only format I could find that gave the actual start positions (and lengths) of each of the pieces that formed a hit (that is, each gap forces a new block).
The regular BLAST output is human readable and is a nightmare to parse. The XML format mentioned has a structure for each HSP, but an HSP can still include gaps:
  [...]
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>51.6</Hsp_bit-score>
      <Hsp_score>122</Hsp_score>
      <Hsp_evalue>9.8e-06</Hsp_evalue>
      <Hsp_query-from>2</Hsp_query-from>
      <Hsp_query-to>98</Hsp_query-to>
      <Hsp_hit-from>35</Hsp_hit-from>
      <Hsp_hit-to>65</Hsp_hit-to>
      <Hsp_query-frame>2</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>24</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>1</Hsp_gaps>
      <Hsp_align-len>32</Hsp_align-len>
         <Hsp_qseq>ITAIGAGLQGPAGCEVIDAGGLLVMPGGIDTH</Hsp_qseq>
         <Hsp_hseq>IAAVGTGLE-PAGAEIIDAGGLLVMPGGIDVH</Hsp_hseq>
      <Hsp_midline>I A+G GL+ PAG E+IDAGGLLVMPGGID H</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

I'm trying to be lazy, I don't want to have to parse the Hsp_qseq and Hsp_hseq to determine exactly where my gaps are. I just want co-ordinates of where my read hits a protein. The PSL format (below) tells me this (tStarts):
match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block   blockSizes      qStarts  tStarts
        match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
---------------------------------------------------------------------------------------------------------------------------------------------------------------
[...]
21      4       0       0       1       8       1       24      ++      <PROTEIN>   662     321     354     <READ>      101     2       101     2       8,17,   321,337,        2,50,

As an alternative, I've executed the BLAT aligner as a means to find alignments of translated DNA reads against a small protein database. I selected BLAT as it provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing). That is, the PSL has a comma delimited field that enumerates each of the starting positions of a query to a target.
However, BLAT outputs a significantly greater number of alignments than alternative tools. This is likely due to its identity percentage threshold defaulting to 25% when emulating a blastx style query.
I can cope with this (apologies to my disks for the I/O), but for filtering purposes, the output PSL format appears to leave me nothing but raw match counts to work with.
The PSL gives me what I want, except now I don't have information on its significance. I know that the BLAT FAQ describes a Perl script for adding "web-based identity and score calculations". I can actually get identity percentages and scores from this script, but I cannot tell what the score actually represents. Is it comparable to bit-score? How can I recover something like an e-value from a PSL output file?

Just out of curiosity: can you elaborate on how BLAT 'provided the easiest means to access information on the actual alignment blocks (with the smallest amount of parsing)'? What about blast HSPs?

@holmrenser Don't you only get the HSPs by parsing the XML output? This is some laze on my part, as the HSPs in my XML described the hits from the perspective of the translated DNA, but I want to know where the alignments on the original sequence are without too much of my own error-prone interference. BLAT just gives the offset. It might be that I've overlooked something in the HSP data though...

No! Not at all! The default output format for blast is a simple text file and there are flags to get one HSP per line (no sequence, only the scores etc). You can get XML as well but that's just a way of making your life harder IMO. Please [edit] your question and clarify what the original (scientific) question you are asking is. Are you looking for exact matches? Distant homologs? Finally, it would be better if you could include an example of your PSL file if that's really what you want parsed.

@terdon Despite my dislike of XML parsing, it seems a more structured way to get out the components of a hit compared to trying to infer structure from a sequence of lines. For what it's worth BLAT doesn't have flags for outputting an HSP per line.

BLAT might not, I admit I've only ever used BLAT online for simple things. BLAST, however (both NCBI and wublast), certainly does. Look at the "tabular" option for the -outfm flag in NCBI blast for example. I've spent more time than I care to remember parsing blast output, if you clarify the exact question you have, I should be able to help.

@terdon Thanks, I know what options BLAST has, but I'm asking about something else. I have updated my question to try and clarify though!

Fair enough (and sorry, didn't mean to appear pushy). I was just wondering if using another tool (like BLAST which would give you the e-values in a parseable table format) wouldn't be a simpler approach. I realize that "redo the analysis with this tool instead" is not a perfect solution though :).

@terdon Hah, no worries. A lot of bioinformatics is "damn, I need to run this with another tool". I'd happily move this to chat if you have any ideas. For what it's worth, BLAST is too slow (for the size of the data) and still doesn't give the granularity on gap information that I'm looking for :( Any thoughts?

It would appear that BLAT and exonerate (as suggested by @terdon) are the only pieces of software that provide enough information for "HSPFragments" to be parsed from their outputs, at least according to BioPython's documentation, anyway:

Most search programs only have HSPs with one HSPFragment in them,
  making these two objects inseparable. However, there are programs
  (e.g. BLAT and Exonerate) which may have more than one HSPFragment
  objects in any given HSP. If you are not using these programs, you can
  safely consider HSP and HSPFragment as a single union.

It would seem the options are sticking with BLAT and generating scores afterwards with the script described in their FAQ, or using exonerate which provides a raw score as part of its default output. What these scores actually mean I will have to investigate, but for now I think the case is closed!

I would suggest using blastx and parsing the output XML. As an example, have a look at this: http://www.perlmonks.org/?node_id=1006367.
The relevant part of the hit would be the following:
<Hit>
  <Hit_num>1</Hit_num>
  <Hit_id>gnl|BL_ORD_ID|1515029</Hit_id>
  <Hit_def>43989.cce_0262 (Cyanothece ATCC 51142)</Hit_def>
  <Hit_accession>1515029</Hit_accession>
  <Hit_len>65</Hit_len>
  <Hit_hsps>
    <Hsp>
      <Hsp_num>1</Hsp_num>
      <Hsp_bit-score>40.0466</Hsp_bit-score>
      <Hsp_score>92</Hsp_score>
      <Hsp_evalue>0.00664016</Hsp_evalue>
      <Hsp_query-from>155</Hsp_query-from>
      <Hsp_query-to>253</Hsp_query-to>
      <Hsp_hit-from>12</Hsp_hit-from>
      <Hsp_hit-to>44</Hsp_hit-to>
      <Hsp_query-frame>-1</Hsp_query-frame>
      <Hsp_hit-frame>0</Hsp_hit-frame>
      <Hsp_identity>17</Hsp_identity>
      <Hsp_positive>27</Hsp_positive>
      <Hsp_gaps>0</Hsp_gaps>
      <Hsp_align-len>33</Hsp_align-len>
      <Hsp_qseq>LRGAICSMEHIEEALGKLKDWARKLIELLLGPR</Hsp_qseq>
      <Hsp_hseq>ITGAVCLMDYLEKVLEKLRELAQKLIETLLGPQ</Hsp_hseq>
      <Hsp_midline>+ GA+C M+++E+ L KL++ A+KLIE LLGP+</Hsp_midline>
    </Hsp>
  </Hit_hsps>
</Hit>

You should be able to get the position information from the following tags: <Hsp_query-from>,<Hsp_query-to>,<Hsp_hit-from>,<Hsp_hit-to>.
To translate back to you original DNA sequence I can imagine you need the <Hsp_query-frame> tag.
And then of course you get the evalue for free: <Hsp_evalue>

An HSP is still allowed to contain gaps, though. So I'd still have to parse one of the qseq or hseq :(

Correct me if I'm wrong, but the tags I mentioned (Hsp_query-from, Hsp_query-to, Hsp_hit-from, Hsp_hit-to) are coordinates on the original sequence.

You are correct (I believe)! But the HSP can contain gaps, so I can't just take those co-ordinates as given to extract the DNA from the read between those positions (as I don't want the gapped nucleotides).

The fact that the read aligns to that position, and that the coordinates are on the original sequence means you actually can. Whether or not you want to do a gapped alignment is a different matter I would think.

Not for what I need - long story short is I'm trying to "flatten" the reads into a pseudo-consensus. I need to know exactly where (what actual position) the gaps are, so I can handle them appropriately. I can't just use the starts and ends of the read as given, because taking that subsequence will include one (or more) nucleotides that should be removed to accommodate the gap(s).


text:
How can I work out gram staining information from the species name?
I have a very long list of prokaryote organism names. I would like to know their gram staining (gram positive or gram negative) information.
Please let me know if there is any database that keeps this information.
I am also trying to find a program that can find gram staining information via organism name only as my organism list is long.
Edit: for each species the NCBI taxonomy lineage is also available

You only have their names, or you also have their genomes?

I have only their names, I can able find their taxonomy lineage using NCBI taxonomy database if that helps me to find the gram staining.

The GOLD database can be used to retrieve a list of species names according to their Gram staining: go to the "Advanced search" page and select the either "gram-" or "gram+" under "Organism Fields" -> "Gram Stain". Then on the results page click on the number under "Organisms". On the results page you can select to have 1000 results per page, meaning that you would need to do some manual work.
I should add that as far as I know this information is based on what people submit when they register their genome projects, so I wouldn't expect this list to be curated or 100% accurate.
If you have access to the genome of the species you are interested in (you could download a representative strain for each species using ncbi-genome-download), you could use traitar, which can predict the gram staining, among other interesting broad phenotypes.

Thanks for your help. I was thinking about downloading GOLD database data and use a program to find my required information. But I am unable to download the data. So can you please help me to do that.

Is there any way to query the GOLD database programatically?

No, their usage policy is in fact quite strict: "Downloading GOLD metadata by any means including web searches, downloadable reports, screen scraping, etc., is strictly prohibited without prior permission, if the purpose is to dissect, alter and store the GOLD data in a database (relational or text) and distribute to third parties/general public/collaborators."

AMR R package can assign gram staining based on taxonomy. See mo_gramstain(x, language = get_AMR_locale(), ...).


text:
What does a question mark ("?") mean in Picard metrics files when I expect a number (integer, float, etc.)
Does anyone know what conditions must be true for Picard to put a ? where I expect a number?
Here is an example from the output of the tool CollectWgsMetrics:
❯ cat sample.wgs_metrics.txt | head -n8 | tail -n2 | verticalize | grep FOLD
FOLD_80_BASE_PENALTY    1.33476
FOLD_90_BASE_PENALTY    2.288159
FOLD_95_BASE_PENALTY    ?

Code searching (link) didn't reveal anything, and googling isn't returning anything obvious either.

[Verticalize is a little CLI tool I've forever been using by Pierre Lindenbaum: verticalize. Easy to compile and I toss it in my ~/bin every time I get a new computer.]

Picard (actually HtsJdk) uses a formatting class (via the MetricsFile class) that emits "?" for NaNs (and Infs).
In your case, there seems to be a division-by-zero error during the calculation of the metric (likely due to zero coverage at some small percentage of bases in your input). This resulting Nan is emitted as a "?".


text:
Multiple sequence alignment with R
I want to incorporate Multiple sequence alignment into our Shiny app, i.e. allow the possibility of selecting different nucleotide sequences, and do a multiple sequence alignment in real-time. The user would define as input different nt. sequences and then a tree will pop up in the app. Do you think this is feasible?
Any recommendations and feedback are really appreciated :)

What do you mean by “tree”? Trees don’t really make sense for pairwise alignments.

I agree with Konrad Rudolph, trees don't make sense with pairwise alignments. There is a pairwiseAlignment function in the Biostring package though, the writePairwiseAlignments funtion can show the alignment.

I mixed up terms here, I meant multiple sequences alignment not pairwise. I have edited the question, sorry about that!

@b.nota Pairwise alignments could be used to compute distances, and a tree could be made from the distance matrix.

@bli Yes, but does it make sense?

@b.nota I think that in theory, a multiple alignment is more apt to capture relevant homologies. So if the tree is meant to represent the phylogeny of the sequences, multiple alignment is theoretically better. That said, there might be practical reasons for which one would have to work with pairwise distances instead.

@bli, okay thanks for your thoughts. I do still agree, however, with Konrad that trees for pairwise alignments don't really make sense.

@SergioArredondo, as the question is closed I cannot provide an answer, but a relevant newer question provides a MWE that shows an outline for the situations where the nucleotides are stored in an msf format, but that can be changed accordingly as other formats can be specified https://bioinformatics.stackexchange.com/a/14160/9311 (and real time is tricky given some of the datasets can be large)

The ‹msa› package on Bioconductor does exactly that.
It doesn’t hand the results on a silver platter, though: you’ll need to read the vignette carefully to learn how to use it. But after that it’s pretty powerful.

"It doesn’t hand the results on a silver platter", that is very true, but some boiler plate code would go a long way since many request the same functionality from a same set of functions


text:
How can I convert SCF trace files to ABI files?
I have a large collection of sequence trace files which are in SCF file format (binary files that in addition to the string of DNA bases also contain the electropherogram of the sample and quality information about the base calls).
I have an existing BioPython pipeline that accepts files in AB1 format (pdf), so I need to convert my SCF files to AB1 files.
I found a question on Biostars that asks for the opposite: converting from SCF to AB1.
The accepted answer there was to use convert_trace from the Staden package.
./convert_trace -out_format scf < trace.ab1 > trace.scf

So I tried the opposite:
./convert_trace -out_format abi < trace.scf > trace.ab1

This did indeed produce a file with *.ab1 extension, but upon importing it to my BioPython pipeline, I got the error
File should start ABIF, not '\xaeZTR'

So apparently, convert_trace converted my files to ZTR format, even altough I specified -out_format abi.
Closer reading of the manual states that [convert_exe] can write CTF, EXP, PLN, SCF and ZTR formats. (so indeed no ABI), but it confuses me why ABI is still a valid output format.
So, my question is, how do I convert SCF files into ABI files?

Apparently convert_trace does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format.
I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by converting AB1 to anything else, often SCF (ZTR came too late and never got enough traction).
That being said, I dug out an SCF I had from 1997, loaded it into Geneious, and re-exported it as .ab1. I was a bit astonished to see that this worked. So Geneious (or similar packages from other companies) may be your best bet for a quick, one-off thing.
Depending on your needs, you might want to think about other ways: either expand BioPython to read SCFs, or maybe to expand the Staden io_lib to be able to write .ab1 (James still maintains that package and is extremely helpful).
Good luck.

I downloaded Geneious and requested a trial license. I'm indeed able to import SCF files, but I don't see 'ab1' as an export format? So, I'm guessing that I will need to go with one of your other two options...

I checked again and my version of Geneious (10.1.3) definitely offers an "ABI chromatogram (*.ab1)" as type in the "File->Export selcted document ..." dialog. Was your SCF selected?

DNA Baser Assembler will convert ABI to SCF (and other formats).
Just double click the ABI file you want to convert then go to the "Samples" menu and choose "Save current samples as...". Plus, it can automatically trim the low quality ends.

There is another tool in that package what will convert all ABI files at batch: Chromatogram Explorer. There is a Lite version that is free, but the conversion to SCF might only work for the first 100 files (if I remember correctly).
Related to the comment posted by BioGeek about the Genious license: DNA Baser Assembler does not require a license key. It will work right away.

The OP was asking for SCF->ABI I think your response is only SCF->ABI.  Do those tools work in both directions?


text:
How to perform functional analysis on a gene list in R?
From an RNA-seq experiment I have about 17000 gene ids for 2 sample conditions arranged according to their log2 fold changes when compared to a control. I need to annotate these, but I've never done annotation before and am wondering how to do this in R.
I'm interested in up and down regulation of the pathways of the involved genes. I'd love to get enrichment and p-value analysis too. I'm interested in whether or not a package similar to PANTHER exists in R for this analysis. I've been reading how different annotation programs can be better or worse; it seems PANTHER is much better than DAVID, so I was wondering if there is a package that provides PANTHER-like analyses in R.
There seems to be multiple packages available, are there any that stand out as being the best?
I'm primarily interested in human samples and annotated pathways.

I've updated your question with that information. What gene IDs do you currently have (e.g., Ensembl, refseq, etc.)?

How do you define best? What do you want to know, the genes that are coordinated, the pathways that are down-regulated or up-regulated? You want a GSEA single-sample or whole sample, competitive or non-competitive? For me this question is like I have a list of number I want to do some operation to get the most important ones, which method is available?

@J0HN_TIT0R to answer a commenter use @ nickname. Please [edit] the question to include this new information.

As it stands, this question is not answerable. @J0HN_TIT0R you don't provide enough background information as to which function you are interested in and why.: this very much decides which enrichment and enrichment tool we can recommend. Also, questions about "the best tool" are discouraged on SE because they are primarily opinion-based.

N.B., I'm avoiding discussion of "best", since that's more or less impossible to answer.
Your question can actually be divided into two:

What's a good tool for pathways analysis (ideally in R)?
What are good sources of pathway information?

For (1), there are a few different possibilities, but I prefer either roast or camera from limma or goseq (a stand-alone bioconductor package). Of those, I expect what roast or goseq are doing are the most similar to what panther is doing in the sense that they're not competitive. However, since genes do correlate with each other, I prefer competitive tests like that offered by camera.
For (2), your options are basically KEGG or reactome. KEGG is problematic, since you need a license for anything remotely recent (there are some R packages to get around this, but I haven't a clue how kosher they are legally). Given that, reactome might be your best bet.
Having written all of that, the best database out there is what IPA uses. This is a commercial product (I have no affiliation), but you can get a demo license that will work for a couple weeks at least. If you're doing analyses relatively often then it makes sense to spend a bit of cash on a license (maybe share with neighboring groups).
Regardless of what you use, you might also want to use pathview, which can produce easily interpretable plots related to pathway enrichment.

I kindly disagree with IPA being best (I thought you wanted to avoid 'best' in your answer). IPA is commercial and is not transparent, but a black box. You have to trust that the company did a good job making this database behind scenes, which is not really how science should be. Just my thoughts.

You requested a tool similar to PANTHER but in R. First, the PANTHER (Protein ANalysis THrough Evolutionary Relationships) tool does a classification  based on evolutionarily related proteins, gene ontologies ( molecular function, and biological process) and pathways.
AFAIK there isn't a tool in R that integrates all these into one, but several packages do pieces of it. However, my observations after using it a couple of times is that is based a lot in gene ontologies. Gene Ontologies (GO) are not pathways, but the best tool I found for GO are topGO (A bit hard to use but useful to get insights on the biology) which test over representations of a term in a given list.

I would go for the topGO package. The logical steps in your pipeline would be:
1) Build a matrix with the following variables:
|gene_ID|GO_ID|GO_term_description|
You can enter a list of uniprot IDs to uniprotkb and retrieve the GO annotations for each gene in your list/file
2) filter that list for the genes that were captured in your differential expression analysis even tho they are not differentially expressed. This filter is defined by DESeq/Sleuth by "at least 4 read counts in at least 49% of the libraries"
3)Define this gene list as gene_universe
4) define your gene set and get the gene_IDs and if possible, also a ranked score (p-val/FDR) of the genes you are interested in
5) plot the graphs and get tables of enriched GO terms (these would be ranked by p-value without need of performing correction analysis
Note: you are going to need some data manipulation functions; I would read the "tidyverse" documentation


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:
Expected allele frequency distribution of SNVs in real NGS data
I have a huge amount of ~20x human WGS samples, aligned, and all SNVs that were called with GATK under standard germline parameters set.
What I need to do is to model SNVs Allele Frequency (AF) for different underlying Copy Numbers. I'd better provide a toy example. For particular genomic region X:
If X is presented by 2 copies for the particular samples, we expect AF to be super-close to 1 or to 0.5.
If X is presented by 4 copies, I expect any particular AF to be close to 0.25, 0.5, 0.75 or 1.
Of course, I can use Binomial Distribution for these purposes. However, as we know, the distribution is not exactly Binomial due to alignment/sequencing biases and the median AF for all heterozygous SNVs is more close to 0.48 but not to 0.5 as we would expect. Another thing: for high copy numbers we expect higher coverages. And GATK use several filters so I suppose that we will not see SNVs with AF like 0.125 (in case if the segment has ploidy 8) - despite the super high coverage there GATK may reject this "weird" AF.
I have read several papers that model SNVs AFs (and I agree that Beta Binomial Distribution may be quite accurate), however, I was not convinced enough that I should use the particular modelling. From your experience (in case if you do SNVs calling), which probabilistic distribution should I use? How should I estimate parameters for each of them (should I expect for CN4 AF=0.5 more frequent than AF=0.75 or vice versa, how to estimate this from data)?
UPD: For simplicity we can say that we have a lot of previously identified regions with ploidy different from CN2, and I can take these coordinates from here. So I can use more or less "supervised" learning for parameters' estimation.

I don't have enough experience to answer which probabilistic distribution should be used.
However, this questions also also asks how to estimate parameters of the distributions. If a binomial distribution is chosen, then Heng Li's paper titled "A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data" 1 is probably the definitive one. Section 2.3.1 of that paper describes an EM algorithm for estimating allele frequencies from multiple samples under the assumption of Hardy-Weinberg equilibrium for arbitrary but constant ploidy.

Ideally even most popular tools like GATK also uses EM algorithm for the record to estimate the likelihood while samtools uses both EM and Brent's method. The real problems with the setting up of ploidy apriori but am not aware as of now about tools that estimates ploidy and then uses that information for SNV calls and finds out the AF. I know tools like ABSOLUTE have the power to estimate the ploidy but then you have to use them as an input for your later SNP calls. But I reckon most tools usually works with the assumption with ploidy set apriori. This is what is coming to my mind as of now.


text:
How to retrieve logical expressions (KO based) for reactions from KEGG?
The completeness of a module can easily be checked by looking at the Definition entry associated with the module. For example, in module M00010, it is given as:
Definition  K01647 (K01681,K01682) (K00031,K00030)

which can be translated to:
K01647 AND (K01681 OR K01682) AND (K00031 OR K00030)

If this expression evaluates to TRUE, the module is considered complete.
I am wondering whether analogous information exists for a single reaction. For example, when looking at reaction R00352, one finds the following information about its Orthology:

But I can't easily find the logical relation that correspond to the KOs for this reaction. It could be:
K01648 AND K15230 AND K15231

or:
K01648 OR K15230 OR K15231

or:
K01648 OR (K15230 AND K15231)

and so on.
In the above example, I know that the correct expression is:
K01648 OR (K15230 AND K15231)

i.e. one either needs either K01648 or both of the other two subunits.
Can this information be retrieved from KEGG for each reaction and if so, how?

Have you tried to asked KEGG directly via http://www.kegg.jp/feedback/ ? If anything, that might give them the information that there's something missing which people need.

You might need to look into the FTP dump; I looked into the KEGG REST API (http://www.kegg.jp/kegg/rest/keggapi.html) and couldn't find anything related to your question...

Some information seems to be available on the module level utilising the API.
Reaction R00352 points to the module M00173. While some information is in the Definition, it's not trivial to link it to the particular reaction. This linkage is however easier accessible via API:

Thanks. There are two issues with this as far as I can see: i) the third KO (K01648) is missing in the list and ii) it is not clear which logical relationship exists between these KOs; it seems it is always an AND but it is not entirely clear (e.g. from the link you provide: K00239,K00240,K00241,K00242,K00244,K00245,K00246,K00247 - do I need all of these or not?.

This table seems to point where the indicated members should all be required as one. In a couple of cases I looked it is indeed the case. Hence ad i) K01648 is separate so should not be noted ad ii) a very detailed decomposition to subunits, but yes.


text:
Determination of intra-run and inter-run repeatability
This question has also been asked on Biostars
Note that this question is specifically in relation to clinical validation and verification activities.
The field lacks established published standards, apart from CLSI documents that are not open-source and cost a non-insignificant sum. However, the technical reproducibility of Illumina sequencing has been well-documented in the literature for over a decade. This consistency is robust enough that this type of analysis primarily measures external factors such as technician differences or day-of-the-week effects. For non-noisy regions, the germline variant-calling congruence is expected to be 100% unless there are significant issues. While this is well-understood, the challenge lies in determining the sufficient level of statistical "modeling" required to demonstrate this rigorously (i.e. statistically sound manner).
Let’s say I have VCF files with germline variants from samples A and B.
Run1: A_rep1, A_rep2, A_rep3, B_rep1, B_rep2, B_rep3,
Run2: A_rep4, A_rep5, A_rep6, B_rep4, B_rep5, B_rep6
I don’t have any “truth” set let us just assume these are random whole-genome patient non-diseased samples.
For analytical performance I want to determine the variability within the run and between the runs per sample.
Does it make sense for me to compare the full pairwise-matrix table and then average it out for precision/recall etc...
Inter-run example (Run1):
Comparisons
A_rep1 > A_rep2
A_rep1 > A_rep3
A_rep2 > A_rep3
But also, I complete the table
A_rep2 > A_rep1
A_rep3 > A_rep1
A_rep3 > A_rep2
Essentially, since we don’t know which false positives (FP) are real, we need to treat them all equally. Ultimately, this will average out the signal, and then the average precision will equal the average recall.
But I don’t really know if this makes sense or how else to approach the problem given no ground truth.
Other alternatives:

Take all shared positives across all 3-samples as “the truth” and
then use that for each sample. But one bad sample will throw all the
rest off (but maybe that’s the point of the measurement). But maybe
this conservative approach is more appropriate?
Jaccard-index

Any thoughts/suggestions from people experienced in this would be greatly appreciated.
Is this the best I can do?
Many thanks

This question needs more clarity; it's hard to tell what you're asking. Please add additional context around what you're trying to do. In particular, make "this question" explicit.


text:
How should I deal with segmental duplications when aligning NGS reads to a reference genome?
This is a follow-up of my other question. I have been having trouble calling variants in the human SMN1 and SMN2 genes, because the human genome has a large segmental duplication there and these two genes (and their surrounding sequences) are essentially identical :

Click here to see the original image on the UCSC browser.
Since variants in the SMN1 gene have been linked to disease, it is important for me to be able to call them in my NGS data. However, because of this segmental duplication, aligning reads to this region correctly is very hard since the presence of the duplication makes each read align almost perfectly to separate locations in the genome, leading to low alignment scores. The low scores will, in turn, block variant callers from calling in those regions.
I thought that the only real solution for this would be to mask the SMN2 gene and its associated duplication in my reference genome which would allow the aligner to align all reads to the region of SMN1. I would then miss any real variants in SMN2, but since the primary gene of interest for  spinal muscular atrophy is SMN1 and not SMN2, I would be OK with that.
My questions are:

Is this a reasonable approach? I was told that "generally, masking segdups is not recommended" by someone who most certainly knows what he's talking about. So is this not a good approach for my stated goal?

If not, then what else can I do that will let me call variants in the SMN1 gene?

I would like something similar to the approach described here: mask one of the regions to rescue mutations in the other.

Short answer for short reads: you can't. Not reliably, at least.
There's a known error mode with aligning short reads to reference genomes where it's quite good at mis-calling reads within repetitive regions. See, for example, this paper, which compares Illumina vs Nanopore for detecting bacterial strains.
If you know that a region is likely to contain a segmental duplication and/or deletion (e.g. due to unexpected changes in read coverage), then it'd be a good idea to double-check any discovered variants (or unexpected non-variants) using a different technology that does not have the same coverage issues. I recommend a long-read sequencing technology that can handle regions at least as large as the segmental duplication of interest (e.g. PacBio up to ~15kb, Nanopore for anything longer), so that you have a high chance of mapping reads uniquely.

I agree that a longer read technique would help but unfortunately that isn't an option for me as I am processing the data generated by others. I need to improve the results of using short reads and don't have the option of using long. But I don't see why you say "you can't". It turns out I most certainly can, and the approach described in the paper you linked to in the question (thanks, I should have included that) works. The complication is dealing with different repeat numbers and ploidy, but the general approach certainly seems to work so far and is giving more accurate results.


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Thank you! I will post a new question.


text:
Map domain names from UniProt bed files to domain accessions
I want to get a bed file mapping human protein domains to the human genome. UniProt actually offers such a thing here. The problem, however, is that the file doesn't include any kind of domain accession, so I have no way of knowing exactly what domain this is. For example:
$ grep P30443 UP000005640_9606_domain.bed | sed -E 's/\t+/ /g'
chr6 29933668 29933929 P30443 0 + 29933668 29933929 153,153,255 1 261 0 . P209-T295; Ig-like C1-type
chr6 29936281 29936542 P30443 0 + 29936281 29936542 153,153,255 1 261 0 . P209-T295; Ig-like C1-type

This tells me there's a "Ig-like C1-type" domain at positions 209-295 of the protein. But what is this domain?
The UniProt page for this protein doesn't add any more information:

I thought they might be using the official description of the domain from PFam. However, there is no such domain in PFam as far as I can tell. I downloaded the full PFam annotation of SwissProt from PFam's FTP site and searched for P30443:
$ zgrep -A4 P30443 swisspfam.gz | sed -E 's/  */ /g; s/[-=_]+//g'
>1A01HUMAN || P30443.1 365 a.a.
C1set 1  (5782) PF07654.14 Immunoglobulin C1-set domain 210290
MHCI 1  (1847) PF00129.17 Class I Histocompatibility antigen, domains alpha 1 and 2 25203
MHCIC 1  (439) PF06623.10 MHCI Cterminus 337364

This tells me there is a "Immunoglobulin C1-set domain" at 210-290. So both the coordinates and the domain "name" are different.
I also checked the Pfam-A.regions.uniprot.tsv.gz file available here. This one had:
$ zgrep P30443 Pfam-A.regions.uniprot.tsv.gz | sed -E 's/\t+/ /g'
P30443 1 8667AFF3F06C4932 d5f7ef1d93bd6e2e0ca1f7b3477ffb67 PF07654 210 290

This gives me a PFam accession (PF07654) but, once more, I can't use this to go back to the bed file because the coordinates are different.
I hoped I could figure this out using the SwissProt flat files, but no luck there either. These are the relevant lines from the entry for P30443:
ID   1A01_HUMAN              Reviewed;         365 AA.
AC   P30443; O77964; O78171; Q9MYA3; Q9TP25; Q9TQP5;
[. . .]
DR   InterPro; IPR021157; Cyt_c1_TM_anchor_C.
DR   InterPro; IPR007110; Ig-like_dom.
DR   InterPro; IPR036179; Ig-like_dom_sf.
DR   InterPro; IPR013783; Ig-like_fold.
DR   InterPro; IPR003006; Ig/MHC_CS.
DR   InterPro; IPR003597; Ig_C1-set.
DR   InterPro; IPR011161; MHC_I-like_Ag-recog.
DR   InterPro; IPR037055; MHC_I-like_Ag-recog_sf.
DR   InterPro; IPR011162; MHC_I/II-like_Ag-recog.
DR   InterPro; IPR001039; MHC_I_a_a1/a2.
DR   InterPro; IPR010579; MHC_I_a_C.
DR   Pfam; PF07654; C1-set; 1.
DR   Pfam; PF00129; MHC_I; 1.
DR   Pfam; PF06623; MHC_I_C; 1.
DR   PRINTS; PR01638; MHCCLASSI.
DR   SMART; SM00407; IGc1; 1.
DR   SUPFAM; SSF48726; SSF48726; 1.
DR   SUPFAM; SSF54452; SSF54452; 1.
DR   PROSITE; PS50835; IG_LIKE; 1.
DR   PROSITE; PS00290; IG_MHC; 1.
FT   DOMAIN      209    295       Ig-like C1-type.

So, while it does have this domain, uses the same "name" and coordinates I had in the bed file ("Ig-like C1-type"), and it has the relevant PFam accession (PF07654), there's no way to link the FT line with "Ig-like C1-type" to the relevant PFam accession.
I even tried using ProSite and downloaded prosite.dat from here. This one has yet another "name" for the domain:
AC   PS50835;
DT   01-MAY-2003 CREATED; 10-MAY-2017 DATA UPDATE; 20-DEC-2017 INFO UPDATE.
DE   Ig-like domain profile.

Finally, I had hoped I could use the coordinates given in the bed file to map but since those aren't the same as those I find in the various PFam files, I'm stumped.
Is there any way to map the domain names given in the UniProt bed files to specific domains (either PFam, or PROSITE or anything else) so that I can know exactly what domain is annotated there?
Update: After jumping through various hoops using the files mentioned above, I have managed to map only 701 domain names which leaves me with 2764 unmapped.


text:
What is the correct way of setting the genotype after splitting multi-allelic sites in a VCF file?
This question was also asked on Biostars
Consider this minimal VCF file which contains a multi-allelic call made by a somatic variant caller on a human tumor sample:
##fileformat=VCFv4.3
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##FILTER=<ID=clustered_events,Description="Clustered events observed in the tumor">
##FILTER=<ID=multiallelic,Description="Site filtered because too many alt alleles pass tumor LOD">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  foo
chr10   43613843        .       GC      TC,CC,TT        .       clustered_events;multiallelic   DP=5180 GT:AD:AF:DP:F1R2:F2R1:SB        0/1/2/3:0,5154,11,15:0.994,0.00229,0.00306:5180:0,2507,4,8:0,2576,6,5:0,0,2585,2595

So, we have a site with three ALT alleles and a genotype of 0/1/2/3.
The VCF specification has this to say about the genotype field:

GT (String):  Genotype, encoded as allele values separated by either of/or|.  The allele values are 0 for the reference allele (what is in the REF field),  1 for the first allele listed in ALT, 2 for the second allele list in ALT and so on.

As I understand the VCF specification, given the ALT field TC,CC,TT, the genotype of 0/1/2/3 is telling me that the reference sequence is absent but there are heterozygous calls for each of the three variants in ALT. This makes sense because this is a tumor sample, so presumably what we have here is three populations of cells, some with the GC=>TC change, some with GC=>CC and some with GC=>TT. Given the very high support for the GC=>TC change, I would assume this is actually a germline variant and there are small populations of tumor cells in my sample with either the GC=>CC or the GC=>TT variant. Alternatively, and quite likely, the latter two variants are sequencing artefacts.
Regardless, I need to split these into one variant per line. As I understand the specification, after splitting, I should be left with three variants and these specific genotypes (I am showing only the genotype here, for simplicity):
chr10   43613843   .   GC    TC   GT 1/0
chr10   43613843   .   GC    CC   GT 1/0
chr10   43613843   .   GC    TT   GT 1/0

Is this actually correct? Unfortunately, standard tools like bcftools produce something different:
$ bcftools norm -m -any foo.vcf 2>/dev/null | grep chr10 | awk '{print $1,$2,$3,$4,$5,$9,$10}'
chr10 43613843 . GC TC GT:AD:AF:DP:F1R2:F2R1:SB 0/1/0/0:0,5154:0.994:5180:0,2507:0,2576:0,0,2585,2595
chr10 43613843 . GC CC GT:AD:AF:DP:F1R2:F2R1:SB 0/0/1/0:0,11:0.00229:5180:0,4:0,6:0,0,2585,2595
chr10 43613843 . GC TT GT:AD:AF:DP:F1R2:F2R1:SB 0/0/0/1:0,15:0.00306:5180:0,8:0,5:0,0,2585,2595

Here, we have sites with a single ALT variant but which have multiple genotypes and that doesn't make sense to me. How can we have a 0/0/0/1 genotype and only one ALT allele? Is bcftools assuming this is a triploid genome?  Is this a bug in bcftools or am I misunderstanding the VCF specification? What is the correct way to display the genotype after splitting this kind of multi-allelic call?


text:
Create heat map that groups genes by expression within cluster in Seurat
I have a set of cells that I am performing Drop-seq on to look at cell expression.  Among my heat maps for gene expression I want to be able to graph them similar to the graph below:

Where the cells are sorted by cluster on the left axis, and have the genes across the bottom.  However I also want to cluster the genes by the expression within the cluster (like how this graph does).  Is that possible in Seurat?
Currently I have created a heat map that looks like this:

It has the clusters listed on the left axis, and small sub-selection of genes on the bottom.  However I want to be able to order those genes by their expression within those clusters, without having to manually order the genes when I insert them into the "feature" argument.
Here is a sample of the code I am currently using to plot the graph:
(Note: I have an object "Combined" that has the non-averaged cell expression, but for now I do want to use average expression)
png("ExpressionHeatmap%03d.png", res=600, width=14, height=7, units="in")
top10 <- Combined.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
avg <- AverageExpression(Combined, features = NULL, add.ident = NULL, return.seurat = TRUE, verbose = TRUE)
Combined.markers %>% group_by(cluster) %>% top_n(2000, avg_logFC) -> top

cc.features <- readLines(con = "/home/alex/lab/drugres/DropSeqFiles/regev_lab_cell_cycle_genes.txt")
s.features <- cc.features[1:43]
g2m.features <- cc.features[44:97]
CCS <- CellCycleScoring(object = Combined, s.features = s.features, g2m.features = g2m.features, set.ident = TRUE)
DoHeatmap(avg, features = c("FTH1P3", "MET", "CYR61", "NFKBIA", "RIPK2", "IFI16", "PLK2", "TNFSF10", "CEBPD", "PNRC1"), group.by = "ident", size = 3, angle = 0, combine = TRUE, draw.lines = FALSE) + theme(axis.text.x = element_text(size = 7), axis.line = element_line(colour = "#ffffff")) + scale_fill_gradient2(low = '#1000ff', mid = "#ffffff", high = '#aa0101', space = "Lab", na.value = "#ffffff", midpoint = 0, guide = "colourbar", aesthetics = "fill") + coord_flip(expand = TRUE, clip = "on")

I tried to use the "FindAllMarkers", "group_by" and "top_n" functions to sort the gene data, and while it helped it still didn't fully sort the genes the way I wanted.
markers <- FindAllMarkers(Combined, features = c(genes-i-am-looking-for), only.pos = TRUE, min.pct = 0.25, thresh.use = 0.25, test.use = "biomod")
geneorder <- markers %>% group_by(cluster) %>% top_n(n = number-of-genes, wt = avg_logFC)

I would then replace the 'features' argument in the DoHeatMap function with the 'geneorder' object.
I know that doing a top10 <- combined.markers style function can get me the effect I'm looking for, but I want to do it, with specifically selected genes, rather than just the top 10 expressed per cluster.
Note: this is for v3.1 of Seurat; the latest version of Seurat at the time of asking.


text:
Acquiring & Processing single cell data (not 10x)
I have a question related to processing of Single cell RNA-seq data originating from non 10x platform. I have worked with data that originated from 10x platform and I parse it through cell ranger pipeline and then put it through Seurat.
As an example the dataset in this GSE156644 repository are not being processed through Cell ranger pipeline and I get an error that indicates "single-end read found" another example is of this GSE115235, in the later I have bunch of files that look like bulk RNA-seq.
If you check the GSE115235 it has bunch of files that are not cellranger compatible. However, they are from single cell experiment. I find it strange that NCBI accepts such submission and the reviewers don't even pay attention to what has been submitted.
I guess, that these are SE RNA seq data that can be aligned using STAR or Bowtie and then merged in some way to get per cell per sample expression.
Along similar lines, how do I deal with 2 fastq files for a scRNA seq run? example https://www.ebi.ac.uk/ena/browser/view/PRJNA602526
This link has scRNA seq files that look paired end (R1 & R2 but NO I1). I am bit confused as cell ranger count option requires the folder to have 3 files (I1,R1,R2). How can I run cell ranger with R1 & R2? Or do I need to use another aligner?
I will send an email to sra and also to the authors, but any other help to navigate through these datasets and process them is deeply appreciated.
Thank you.

Cell ranger does not require I1 and it was a bad idea for them to ever support that. You can directly run cellranger count on demultuplexed fastq files. 2. Salmon alevin or STAR-solo would be examples of common alternatives. Salmon alevin is in fact much faster and likely better.

Modified from @devon-ryan's comments:
The couple of samples I looked at have the sequence from read 1 embedded in their header. I1 isn't needed, they're already demultiplexed. You could feed those into STAR or similar tools, but processing them will be more of a pain since most tools expect the original unprocessed version of what they uploaded.
Modified from atpoint's comments:
Single-cell data is not really useful as "bulk" as the notable amount of PCR during the library prep creates quite a bias that without the UMI information cannot be corrected. I'd email the authors and ask for the R1 files. Single-cell uploads at NCBI tend to be a mess, that's unfortunately how it is. Without R1 data are almost useless.
Those file format issues are strange, but it is what it is. Reviewers do not pay attention to method sections, choice of tools and data submission most of the time in my experience, especially not when the paper is biology-focused rather than a computational-focused one. Even more, you often do not find the suitable metadata for a proper analysis such as information about batches etc. The only choice you have is emailing the authors as noted above.
If by "two files" you mean that the same library (=sample) was sequenced over multiple lanes (very common to do that since required depth is often larger than what a single lane can provide) then merge the fastq files, so all R1s into one and all R2 into one file.


text:
Why does my HWE QQ plot have extreme deviation and what does it mean?
This question was also asked on Reddit
I have recently completed my thesis and one of the comments was that I report on why this graph looks this way. I have tried to find a reason but the closest I can come to is that it means that there are false positives? Or that this is possibly due to population stratification?
The trait of interest was alcohol dependency with covariates of sex and age, in a non african american population.
I am also unsure as to how this will impact my research or the implications of this?

I don't understand why this population (non-African American) "reacted" this way - do I provide a different null model? The previous populations had "relatively normal" plots. They were all run through the same pipeline.
How the data was calculated:

There were 1061 individuals in the non-African American dataset. At QC phase 0, there were 2 individuals who had discordant sex information and 87 duplicate SNPs. When the MAF value was set to 0.000015, 125708 SNPs were removed when the MAF. Both the mind and geno parameters were set to 0.2, resulting in 127357 and 0 SNPs to be removed respectively. When setting the HWE to $1 \times 10^{-12}$, 4495 SNPs were removed. After all these parameters were set and the pipeline completed, there were 851773 SNPs left from 1044 individuals.

This much deviation from the expected distribution looks like a modeling issue, i.e. that the statistical model use to represent the null hypothesis is incorrect.
I can't see any values that lie along the expected line. It's more common to see a generally linear trend, with a few outliers just at the end.
It could very likely be due to population stratification. You have defined this population as "non-African American", which seems to me to be quite a vague definition that would capture a lot of different population groups with distinct genetic histories.
In my own research many years ago, I found substantially different haplotypes within the ADH gene cluster between Māori and European individuals, and I could imagine that not adjusting for population structure in other populations would have an impact on association scores.

The level of discordance is extreme, particularly for humans. This is either a publication in Nature/Science and massive advance in the field (especially humans), or the genetic factors have massive confounders.
In theory, I'm aware there must be a genetic component to alcoholism, there is a biochemical component about the rate of processing alcohol toxins. However, non-genetic (especially social/community) affects of alcoholism could predominate.
What I personally think is happening is there are very strong non-genetic confounders which are disrupting the model and this is getting confused with neutral HWE. For example, if alcoholism predominates in a community or social group the transmission is social not genetic, so the confounder is extreme because they are genetically related, but the gene(s) are not being inherited. In phylogenetic terms it will make the "gene" appear far higher up (far deeper) in the tree than it actually is. The "lone" drinker is probably the model, but group behaviour and peer-pressure is probably more realistic. Equally community effects could be suppressing alcoholism, such as religious belief.
Thus if a non-genetic model of alcoholism, (loads of social factors) resulted in a stronger fit by for example 3 logs difference, whilst the genetic model is 50 logs difference, then that is the answer and points to the need for careful population sampling for GWAS. Upvotes J Human Gen is cool.

As a disclaimer, I do not understand not human genetics outside immunology. I do understand multiple linear regression models and failing to account for a pivotal factor(s) causes unusual effects. An obvious example in this case is religious belief: thus alcoholism genes could be present but the environment is not presented. For example, you might conclude long hair caused a reduction in human growth because they tightly correlated in humans. However, once gender is a factored into the model, the result is clear hair-length associates with gender not height.


text:
Coverage required
I was came across a problem during an exercise in a book and I don't really know how to solve it. I feel like something's missing.
"coverage, c = $NL/G$ (N=number of reads, L=read length, G=genome length)
total expected gap length = $G(e^{-c})$
total number of gaps = $N(e^{-c})$
Question: You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from equimolar base composition?"
The answer from the book is basically just: Solve $N(e^{-c})=4$ for $c$, but the question doesn't give what $N$ is?

What do you mean? Isn't N the number of reads? That's what you say in your question.

yeah, what I meant was that the question doesn't give me the value of N, without which I can't solve for c

Solving for N and solving for C are the same thing (you know L and G in the first equation after all). It also doesn't give you a nice answer because the actual value (11,926,375 reads, which is ~14.9X coverage) isn't really a nice pretty number.

According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113
(a) What fraction of a genome could you expect to assemble from eightfold coverage?
(b) What total gap length would you expect in an assembly of a 2 Mb target genome size from eightfold coverage?
(c) How many gaps would you expect in an assembly of a 2 Mb target genome size from an eightfold coverage of fragments with a read length of 500?
(d) You want to sequence a 4 Mb genome by the shotgun method, by assembling random fragments with read length 500. What coverage would you require, to expect no more than four gaps, assuming no complications arising from repetitive sequences or far-from-equimolar base composition?

I'm agreed that whitout knowing $ N $ it's not obvious to answer this last question so perhaps they're a mistake in the statement of this exercise.
As Devon Ryan suggested, let's assume eightfold coverage means the coverage $ C $ = 8X
So we can calculate the number of reads:
$ N = \frac{GC}{L} $
$ N = \frac{4000000\times{8}}{500} = 64000 $ reads
Then we can calculate the required coverage to expect no more than four gaps:
$ N( e^{-C}) = 4 $
$ e^{C} = \frac{N}{4} $
$ C = log{\frac{N}{4}} $
$ C $ = 10X
(PS: I would like to thank @LordOfRacoons to help answer this question.)


text:
Read alignment using Bowtie2
So this is related to CRISPR-CAS9. I am working with off-target predictions for my thesis and was looking at all scientific papers related to CRISPR. I found one and decided to use their datasets. The datasets were not available directly and the method used to obtain them was specified which I had mentioned below:
They have around 30 sgRNAs which are each 23bp long and are aligning the sgRNAs to human reference genome with a tolerance of upto 6 mismatches using Bowtie2.
In order to access the exact dataset, I thought I need to replicate the method they used. To do so, I collected the sgRNAs from their paper and tried to use Bowtie2 to align them against human genome with upto 6 mismatches.
Since my laptop does not have the capacity to perform this operation (it has 8GB RAM), I am searching for interfaces that provide access to Bowtie2. I found Galaxy to be useful. Does anyone know of any other tool that provides an option for Bowtie2?
Thanks!

One possible workaround would be to blast your sequences against the genome. Then, when you have found where they align, you can extract those regions only into a new fasta file (make sure to take a few kb around each target region) and then use that fasta file as the "genome" you pass to bowtie2.
That would give bowtie a far smaller file to work with and should reduce the memory needed significantly.

This would be a very good alternative. Thanks a lot @terdon for the answer.


text:
ATAC-seq macs2 peak splitting in sliding windows
This question has also been asked on Biostars
I used macs2 to call peaks for atac-seq data. now my goal is to split the peaks into 50 bp windows with 25 bp steps and then calculate the Tn5 integration frequency in each window.
How should I proceed with that?

import numpy as np
np.random.seed(0)
import pyranges as pr

gr = pr.random()
gr.Score = np.random.randint(100, size=len(gr))

gr = gr.slack(25) # make data wider for this example

print(gr)

t1 = gr.tile(50)

def increase_by_25(df):
    df = df.copy()
    df.Start += 25
    df.End += 25
    return df

t2 = t1.apply(increase_by_25)

tiled = pr.concat([t1, t2]).sort()

print(tiled)
# +--------------+-----------+-----------+--------------+-----------+
# | Chromosome   | Start     | End       | Strand       | Score     |
# | (category)   | (int32)   | (int32)   | (category)   | (int64)   |
# |--------------+-----------+-----------+--------------+-----------|
# | chr1         | 5205300   | 5205350   | +            | 33        |
# | chr1         | 5205325   | 5205375   | +            | 33        |
# | chr1         | 5205350   | 5205400   | +            | 33        |
# | chr1         | 5205375   | 5205425   | +            | 33        |
# | ...          | ...       | ...       | ...          | ...       |
# | chrY         | 41326450  | 41326500  | -            | 3         |
# | chrY         | 41326475  | 41326525  | -            | 3         |
# | chrY         | 41326500  | 41326550  | -            | 3         |
# | chrY         | 41326525  | 41326575  | -            | 3         |
# +--------------+-----------+-----------+--------------+-----------+
# Stranded PyRanges object has 7,964 rows and 5 columns from 24 chromosomes.
# For printing, the PyRanges was sorted on Chromosome and Strand.

This is as far as I can get without example data and a clearer explanation of what you need.