text:
Classification (supervised learning) of expression data on pathway level
I was wondering if there is any way to apply classification algorithms (e.g random forest) on microarray data but not using the genes as predictors/features but the pathway they belong to.
The thing is that the expressions of the genes that belong to the same pathway should somehow, grouped together and transformed into a new value that will describe that pathway for each specific sample.
Let me give an example.
Assuming that we have g1, g2, g3, g4, g5, g6, g7, g8, g9, g10 as genes
and p1, p2, p3 as pathways.
And g1, g2, g3 belong to p1 while g3, g4, g5, g6belong to p2 and g2, g6, g7, g8, g9, g10 belong to p3
One thought is to calculate the mean expression of the above gene groups per sample and then run the classification algorithm on top of these values. But I would like to hear/read your opinion on such an approach and also your suggestions on using a different method to do this.
The final goal would be the extraction of the most significant classification predictors (pathways).

Do you assume that expression of genes are similar in pathways? Because in a real experiment they are not.

No, I don't want to assume something like that but I mentioned the mean in order to group the genes of each pathway. This approach might be totally wrong but that's why I posted here. To get feedback and different opinions.

My point is, that I haven't seen an algorithm that can transform the expression values of different genes into one value that represents a given pathway.

Oh ok. I see...

It looks like you're trying to incorrectly reinvent the activation Z-score produced by IPA or similar tools. I suggest you simply use that instead.

IPA is commercial. Isn't it? Can you provide me with free R-based ones to understand what exactly you suggest?

Have a look at the GSVA package. It allows to convert a matrix with genes x Samples to a pathways x Samples using several methods ssgsea, gage, gsva...
Afterwards you can use that matrix as input for differential expression of pathways or classification algorithms or whatever.

However it depends on the input of the "pathways" you give it. Make sure that the pathways are correct. Also consider that the pathway topology might change with conditions.


text:
Why is my STAR reference genome indexing aborting on my GNU/Linux server but not on my Mac OS X laptop?
I am running the following command to index my genome:
STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gff3 --sjdbGTFtagExonParentTranscript Parent --sjdbOverhang 100
and I get the following output:
Feb 27 10:27:37 ..... started STAR run
Feb 27 10:27:37 ... starting to generate Genome files
Feb 27 10:27:40 ... starting to sort Suffix Array. This may take a long time...
Feb 27 10:27:42 ... sorting Suffix Array chunks and saving them to disk...
Feb 27 10:28:26 ... loading chunks from disk, packing SA...
Feb 27 10:28:30 ... finished generating suffix array
Feb 27 10:28:30 ... generating Suffix Array index
terminate called after throwing an instance of 'std::bad_alloc'
  what():  std::bad_alloc

Searches for the GNU/Linux error lead to some pages about C++ that I don't understand.
Submitting it as a batch job works fine but I would like to know why it does not work on a node. It is quite practical for me to test small code like this without submitting a batch job.
The command above works fine on my laptop running Mac OSX.
The STAR version is STAR_2.5.3a on both computers: STAR_2.5.3a
The memory usage of the node is the following:
total            used       free     shared    buffers       cached

251              251         0         0       0             232

so, is there 0 memory for me to use on this node?
Why is my indexing aborting?

Are they the same version of STAR? Compiled with the same flags? Do you have enough RAM on the node you are using?

What do you mean by 'compiled with the same flags'? I have installed STAR using conda from the bioconda channel. Do you know how I can check the available RAM on my node?

That will depend on how your system is set up, but if you have ssh access to your node, you can just run free -h (or just free) to see available memory. That said, you usually don't want to be running things directly from the nodes. Most such systems will have a queuing tool that will distribute jobs across nodes and that will help you avoid this type of problem, and ignoring the queue means i) you may well have the kind of problem you are encountering and ii) you are being rude to the other users of the system; you are effectively cutting the queue.

A std::bad_alloc error tends to mean that you've run out of memory. It's not unusual for head nodes to be fairly limited in capacity, so presumably that's the issue.

source: Stack Exchange

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

sort: votes


text:
Create galaxy local toolshed using Planemo
I tried to create galaxy local toolshed using planemo. I used command as follows,

create virtual env
virtualenv .venv; .venv/bin/activate

Install planemo
pip install planemo

To configuring a shed account
planemo config_init

Create a repository
planemo shed_init --name=seqtk_seq \
                  --owner=planemo \
                  --description=seqtk_seq \
                  --long_description="Tool that converts FASTQ
                                      to FASTA files using seqtk" \
                  --category="Fastq Manipulation"
planemo shed_create --shed_target local

This produces the following error:
Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.
(.venv) varuna@varuna-Inspiron-15-3567:~/git/galaxy_stable/planemo

Hi, could you paste the text of the error instead of screenshots? It might make it easier to read and help you.

Failed to include files for [{'source': '**'}]
Problem encountered executing action for one or more repositories.

Could not update local
HTTPConnectionPool(host='localhost', port=9009): Max retries exceeded with url: /api/repositories?key=dc5c7cd5cf56a2f44aacf755109e63b9 (Caused by NewConnectionError('<urllib3.connection.HTTPConnection object at 0x7f55542b7310>: Failed to establish a new connection: [Errno 111] Connection refused',)), 0 attempts left:

Try providing a project name planemo shed_create --shed_target local someProjectName.

I created galaxy local toolshed editing the xml files inside galaxy folder and  ran it using 9009 port. After run toolshed, I tried above mentioned command with project name.
planemo shed_create --shed_target local /home/varuna/git/galaxy_stable
Got following error.
"description required for automatic creation or update of shed metadata"

--shed_target local will only work if you're running a Galaxy toolshed on your local computer, since it's reading the .shed.yaml file you made with planemo shed_init and trying to create a repository with that name or a local (in this case) server. Planemo itself won't start a Galaxy toolshed, since it's fairly unusual for people to want to do testing involving creating local repositories. What you probably want to do instead is use planemo test instead. Once the planemo testing passes then use the various shed* commands to upload your wrapper to the toolshed or testtoolshed.

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


text:
gffread: GFaSeqGet errors on coordinate overhang
Disclaimer: I had this issue posted on Tuxedo Tools users group and shared it on Twitter, but could not get an answer to this from the developers, nor find documentation of this issue online. So, I'll share my solution here below and see if some of you have a better solution to add.
Using gffread utility to convert transcripts from GTF as output of Cufflinks to FASTA, I get some GFaSeqGet errors, as follows (here is the full log file):
$ gffread transcripts.gtf -g ../genome.fa -w transcripts.fasta
Error (GFaSeqGet): end coordinate (32519) cannot be larger than sequence length 32512
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (7812) cannot be larger than sequence length 7800
Error (GFaSeqGet): end coordinate (6327) cannot be larger than sequence length 6326
Error (GFaSeqGet): end coordinate (4651) cannot be larger than sequence length 4648
Error (GFaSeqGet): subsequence cannot be larger than 3551
Error getting subseq for CUFF.29386.2 (1..3568)!

That's a bit weird since I used the same genome file as the one used to generate the alignment file required by Cufflinks, and the coordinates overhang does not make a lot of sense to me.
The alignment file (SAM) was generated using Hisat2:
$ hisat2
    --rna-strandness RF
    --dta-cufflinks
    -x ./hisat_idx
    -q
    -1 ../1A_1.fq.gz,../1B_1.fq.gz,../2A_1.fq.gz,../2B_1.fq.gz,../3A_1.fq.gz,../3B_1.fq.gz
    -2 ../1A_2.fq.gz,../1B_2.fq.gz,../2A_2.fq.gz,../2B_2.fq.gz,../3A_2.fq.gz,../3B_2.fq.gz
    -S out.sam

Then sorted:
$ cat out.sam | sort -k 3,3 -k 4,4n -T . > out_sorted.sam

My question is: are only the few sequences reported in the log affected? Or is it a more general problem? The answer seem to be the second one.
Out of gffread, alongside the log file, I actually got a FASTA file with no additional error messages during its generation (537M out of 404M GTF file) - but there are around 13k transcripts which are present in the GTF but seem to be missing from the FASTA:
$ cat transcripts.gtf | awk '$3=="transcript"{print $0}' | wc -l
231139
$ cat transcripts.fasta | grep '^>' | wc -l
218075

Do you know what would be a cause of this problem?

Basically, the source of the problem is issues with hard/soft clipped bases in short-reads RNAseq reads alignments to reference genome in post-transcriptome assembly steps.
Truncating the coordinates in the gtf file to match the reference is a solution. For this, I will share a perl solution (non-official) that will do the job.
perl cufftrim.pl genome.fa.fai transcripts_to_fix.gtf > transcripts_fixed.gtf

Where:

genome.fa.fai is the fasta index of the genome file
transcripts_to_fix.gtf is the transcript set out of the assembler with issues
transcripts_fixed.gtf is the new transcript file with the repaired end coordinates

Content of cufftrim.pl:
#!/usr/bin/perl
use strict;
use Getopt::Std;

my $usage = q/Usage:
 cufftrim.pl genome.fa.fai transcripts.gtf > transcripts.trim.gtf
/;
umask 0002;
getopts('o:') || die($usage."\n");
my $outfile=$Getopt::Std::opt_o;
if ($outfile) {
  open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
  select(OUTF);
  }
# --
my $fai=shift(@ARGV) || die ("$usage\nSamtools fasta index required!\n");
my $gtf=shift(@ARGV) || die ("$usage\nCufflinks' transcripts file required!\n");
my %slen;

open(FAI, $fai) || die("Error opening $fai\n");
open(GTF, $gtf) || die("Error opening $gtf\n");
while (<FAI>) {
 chomp;
 next if m/^#/;
 my @t=split(/\t/);
 die("Parsing error: $t[1] length for sequence $t[0] ?\n")
    unless $t[1]>10;
 $slen{$t[0]}=int($t[1]);
}
close(FAI);

# now process the GTF file
my ($tcount, $ecount); #trim counts
while (<GTF>) {
 my @t=split(/\t/);
 my $l=$slen{$t[0]} || die("Error: contig $t[0] not found in the fasta index!\n");
 if ($t[4]>$l) {
   $t[4]=$l;
   $t[2] eq 'exon' ? ++$ecount : ++$tcount;
 }
 print join("\t", @t);
}
close(GTF);
print STDERR "$tcount transcripts, $ecount exons adjusted.\n";
# --
if ($outfile) {
 select(STDOUT);
 close(OUTF);
 }

Filtering the SAM file prior to assembly to exclude clipping might also be an alternative solution.


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

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

@terdon Whole genome 30X on HiseqX. Link

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

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

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

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

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

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

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

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

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


text:
coexnet package for network co-expression analysis, is looking for a filelist.txt that doesn't exist
On my way to create a co-expression matrix for a microarray experiment, using the coexnet package I had to use a function called getAffy() which according to the packages' manual it :

searches in your current or designated path file the folder with the
  GSE accession number and reads the filelist.txt file that contains the
  name of each GSM sample, in order to recognize them and join them in
  an unified AffyBatch object.

So in some way, this file is needed to create the AffyBatch object. The thing here is that there is no such a file downloaded from the getInfo() function.
Googling out the filelist.txt I reached only this result from GEOQuery manual, referring to that file. And as you can see there is a command in it eList2 <- getGEOSuppFiles("GSE11675") that is believed that it returns also the filelist.txt. I tried even with the example of the previously mentioned command, but no filelist.txt file ever fetched.
Is there any way to setup my own filelist.txt file? Does anyone know what it contains and in what form?

Finally, I found it. filelist.txt file is accessible from the FTP server of the NCBI. For example for that specific microarray is stored here
So there might be a problem at software side and they don't get downloaded.

I know you now solved the issue, but I think you should double check the real manual of GEOquery. But I think that the problem is with the type of data that is provided, I recall having downloaded this file with getGEOSuppFile


text:
Identifying outlier loci in R
I have done some genome scan analysis and I want to identify the outlier loci. For an expected value of the FRD equal to q = 10, a list of candidate loci can be obtained by using the Benjamin-Hochberg procedure as follows:
q=0.1
L = length(adjusted.p.values)
W = which(sort(adjusted.p.values) < q * (1:L) / L)
candidates = order (adjusted.p.values)[w]

and this is part of the outliers and my output
head(candidates)
 [1]  3  4  9   10  13

But I need the output in a different type.
I want to assign TRUE if a locus has been chosen as outlier and FALSE if it has not
  head (desired_output)

  [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE

Any idea how can I modify the Code?

As Devon said, removing the which will give you W as a logical vector that can be used to index the sorted p-values. To get a logical vector that can be used to subset the original (unsorted) values you’ll need to go a step further:
candidates = is.element(1 : L, order(adjusted.p.values)[W])

However, I’m not sure this will work with your data: the name adjusted.p.values suggests that these values are already FDR-adjusted. So you cannot/should not perform additional correction on them. If they are already adjusted, then you can get candidate loci simply by
candidates = adjusted.p.values < q

Thanks@Konrad Rudolph.candidates = is.element(1 : L, order(adjusted.p.values)[W]) did exactly what I wanted.

Remove the which(), which is converting the boolean vector to a vector if indices to TRUE.


text:
Can I export the SNPs between whole genome alignments using the command line?
I am using Mauve to align two whole genomes. I use the following command to get the alignment in xmfa format:
/Applications/Mauve.app/Contents/MacOS/progressiveMauve ref1.fa ref2.fa --output alignment.xmfa
I then open the xmfa file in the Mauve GUI and export the SNPs. Is there a way of doing this using a command? I have tried to follow the procedure described here but the link to snapshot releases of Mauve does not work anymore.
I could not find anything in the Mauve command line help:
/Applications/Mauve.app/Contents/MacOS/progressiveMauve --help
nor on the Mauve website.
I tried the procedure described by Konrad Rudolph. I wrote:
MAUVE_DIR=/Applications/Mauve.app/Contents/Java
in the terminal. Is this correct? Don't I need an export and quotation marks to create a bash variable? I tried adding an export and using quotation marks for the file but the error below remained the same.
followed by:
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
followed by:
java org.gel.mauve.analysis.SnpExporter -f alignment.xmfa -o aln.snps
and it seems to work.

It works for me with the download from download page. The download archive contains Mauve.jar so the command line in the answer you link to should work.
However, you might need to include all the jar files (not just Mauve.jar) in the classpath for this to command work, though. For instance, the following works for me:
MAUVE_DIR=‹the-dir-containing-Mauve.jar›
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
java org.gel.mauve.analysis.SnpExporter -f aln.xmfa -o aln.snps

… there’s probably a smarter way to make this work but I don’t use Java enough to know how.

Should I download this file: http://darlinglab.org/mauve/snapshots/2015/2015-02-25/MacOS/Mauve-snapshot_2015-02-25.dmg ? But this is to install mauve, right? Where can I find the Mauve.jar file?

@charlesdarwin macOS applications are regular directories. Open the DMG and go inside the Mauve.app dir. Inside that directory you’ll find /Contents/Java/Mauve.jar; to display the contents in Finder, use the instructions here: https://apple.stackexchange.com/a/95942/2392 — but in the terminal you can simply cd into the directory.


text:
How to use C++ htslib to read VCF contig name and size?
A typical VCF file has:
##contig=<ID=chr1,length=248956422>
##contig=<ID=chr2,length=242193529>

I would like to use htslib in C++ to read it. My attempt:
htsFile *fp = bcf_open("my.vcf", "r");
bcf_hdr_t *hdr = bcf_hdr_read(fp);

In https://github.com/samtools/htslib/blob/develop/htslib/vcf.h, I'm not able to find a function that can do that for me.
How to read chr1 and 248956422 in C++?

As far as I can tell, the (internal) function bcf_hdr_parse_line does that. It returns an object x of type bcf_hrec_t*, for which strcmp(x->key, "contig") == 0) for your lines. You can then access the x->keys and x->vals. Untested.

@KonradRudolph Can I use it externally properly? Without hacking?

Depends on what you mean by “properly”. Since the API isn’t properly documented anyway, the distinction between external and internal functions is somewhat fluid. It might break with future versions of htslib, if that’s what you’re concerned about.

@KonradRudolph Yeah. Breaking with future versions is my concern. If I see don't any other alternative, I might go with that or just parse the file myself.

#include "htslib/vcf.h"

int main(int argc, char *argv[])
{
    htsFile *fp;
    bcf_hdr_t *hdr;
    bcf_idpair_t *ctg;
    int i;
    if (argc == 1) {
        fprintf(stderr, "Usage: print-ctg <in.vcf>\n");
        return 1;
    }
    fp = vcf_open(argv[1], "r");
    hdr = vcf_hdr_read(fp);
    ctg = hdr->id[BCF_DT_CTG];
    for (i = 0; i < hdr->n[BCF_DT_CTG]; ++i)
        printf("%s\t%d\n", ctg[i].key, ctg[i].val->info[0]);
    bcf_hdr_destroy(hdr);
    vcf_close(fp);
    return 0;
}

On stability: this use has been in htslib forever. In general, functions/structs/variables in the public headers are meant to be stable. However, there is no guarantee that future versions will always keep the same APIs.


text:
Converting VCF file to PLINK bed/bim/fam files
I am trying to find the best way to convert VCF files to PLINK binary bed/bim/fam files, but it seems like there are many varied ways to do this. (For example, using Plink 1.9 --vcf tag, bcftools, GATK, and vcftools).
Obviously they all probably have their uses and you could likely spend a while discussing this, but can anyone give a short, simple explanation as to whether the methods are generally the same, or whether one or another is generally best for certain situations? Thanks!

Given that Plink reads in VCF files natively, and Plink is your desired output format, I expect that plink --vcf <file> would be the best option.
As mentioned by @sarah, Plink uses the sample IDs to create family/individual IDs, see  "Variant Call Format" section on this page:

VCF files just contain sample IDs, instead of the distinct family and
within-family IDs tracked by PLINK. We offer three ways to convert
these IDs:
--double-id causes both family and within-family IDs to be set to the sample ID.
--const-fid converts sample IDs to within-family IDs while setting all family IDs to a single value (default '0').
--id-delim causes sample IDs to be parsed as <FID><delimiter><IID>; the default delimiter is '_'. If any sample ID

does not contain exactly one instance of the delimiter, an error is
normally reported; however, if you have simultaneously specified
--double-id or --const-fid, PLINK will fall back on that approach to handle zero-delimiter IDs.
If none of these three flags is present, the loader defaults to
--double-id + --id-delim.

There are caveats associated with the types of variants that Plink will pull in from the VCF files, mostly related to limitations of Plink's own file formats. But if your end-goal is Plink, then that should be fine.

Additional tidbit I learned in the last year that I wanted to share with anyone else that's working on this.
Plink creates an extremely generic fam file for you, but if you are updating this fam file to reflect your population's family structure, sex, etc., note that when using plink binary files (bed/bim/fam) the order of the individuals in the fam file needs to match the order of the individuals in the other files (which is the order of the individuals in the generic fam file)
Do be careful if you care about which alleles are major/minor, as sometimes the order gets switched around in plink. For more info see "Variant Call Format" section on this page.

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


text:
Generating the reconstructed alignment from BAM
I have a (small) BAM file with CIGAR and MD fields.
Question 1: What tools exists in Python and/or R to reconstruct the alignment between the reference and the read in a BAM? Given that this is a very standard analysis, I feel that there should be R packages or Python libraries which have functions which do this task given a BAM input...
My goal was to create a *tsv file with two columns, one with the reconstructed alignment and the other column with the corresponding sequence. This would be useful in my analysis. (I would either use R data.table, or pandas in Python.)
Here is read and reference:
CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
|||||||||||||||| |||||||||||||||||
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG

Here is the desired output dataframe (with only the one example row). We could use an R data.table/data.frame, or a pandas DataFrame. Column1 are the alignments (i.e. the reconstruction of the alignment between the reference and the read based on the CIGAR and MD strings), Column 2 are the corresponding sequences. :
ReadAligned    Reference
CGGGCCGGTCCCCCCC-GCCGGGTCCGCCCCCGGG    CGGGCCGGTCCCCCCCCGCCGGGTCCGCCCCCGGC
...    ....

Question 2:
Using the CIGAR/MD tags, one could then create columns with the number of insertions and the number of deletions. This problem feels much more straightforward, e.g. the python library cigar or the cigar parsing functions of pysam

Could you clarify your question: You want to store the alignments in a matrix? Or the CIGAR format in a matrix?

@Llopis A tsv/csv file with two columns: the first column has the alignments, the second column has the sequence. Does this make sense?

Do you want sam2pairwise to generate the alignments for you, or do you want to parse the cigar strings yourself in R/python? Writing a cigar parser is tedious, but do-able. It might be easier if you just use sam2pairwise to output the alignments for each read, then you parse that into a table.

Now your question is reduced to "how to count the number of insertions and deletions" in cigar, which is much easier. With perl, it is just an one-liner echo 78M2I10M5D30M|perl -ne 's/(\d+)([ID])/$h{$2}+=$1/eg; print "$h{I}\t$h{D}\n"'. Someone will probably post python answers.

@user172818 I've tried to edit this to be more specific. Is this more organized? The second problem is easier; I'm primarily motivated by the first problem, which is related.

If all you want are the number of insertions/deletions, you can use this to extract that info:
from cigar import Cigar
c = c=Cigar('20M5I10M3I5M')
num_insertions = len(list(filter(lambda x: x[1] == 'I', c.items())))
insertion_length = sum([x[0] for x in list(filter(lambda x: x[1] == 'I', c.items()))])

Deletions would be the same thing with 'D' in place of 'I'

Is there a standard method to grab the string from the BAM within python? It looks like this is parsed automatically with pysam

I guess pysam would be the standard way

I am answering question 1. However, as I am unfamiliar with python, I am providing a javascript solution only to show the basic logic. To try the code at the end, you need to install node.js and run it with node this-script.js.

Given that this is a very standard analysis

I tend to believe this is a very rare analysis. I have only seen it done twice. Converting CIGAR+MD to padded alignment is complicated in coding, as you have to keep sequence, cigar and MD in sync all the time. The code below generates two strings in one go. It may be logically easier to generate the reference and the query strings separately.
In general, MD is ill defined. It is too difficult to work with. If possible, avoid this tag. You can do a lot of things with the combination of CIGAR and NM.
function md2pad(CIGAR, MD, SEQ) {
    var m, re_cigar = /(\d+)([MIDSHNX=])/g, re_MD = /(\d+)|(\^[A-Za-z]+)|([A-Za-z])/g;
    var cigar = [];
    while ((m = re_cigar.exec(CIGAR)) != null) // parse CIGAR into an array
        if (m[2] != 'H') // do nothing for hard clipping, as it has no effect
            cigar.push([parseInt(m[1]), m[2]]);
    var k = 0, sx = '', sy = ''; // k: k-th cigar operator; x: ref; y: query
    // cx/cy: start coordinate of cigar[k] on ref/query; mx/my: start of current MD op
    var cx = 0, cy = 0, mx = 0, my = 0;
    while ((m = re_MD.exec(MD)) != null) {
        if (m[2] != null) { // deletion from the reference
            var len = m[2].length - 1; // length of deletion
            sx += m[2].substr(1), sy += Array(len+1).join("-");
            mx += len, cx += len, ++k; // consume a deletion D from ref
        } else { // copy or mismatch
            var ml = m[1] != null? parseInt(m[1]) : 1; // length of this MD operation
            while (k < cigar.length && cigar[k][1] != 'D') {
                var cl = cigar[k][0], op = cigar[k][1];
                if (op == 'M' || op == 'X' || op == '=') {
                    if (my + ml < cy + cl) { // an MD ends in the middle of cigar M
                        if (ml > 0) {
                            if (m[3] != null) sx += m[3], sy += SEQ[my];
                            else sx += SEQ.substr(my, ml), sy += SEQ.substr(my, ml);
                        }
                        mx += ml, my += ml, ml = 0;
                        break;
                    } else {
                        var dl = cy + cl - my;
                        sx += SEQ.substr(my, dl), sy += SEQ.substr(my, dl);
                        cx += cl, cy += cl, ++k; // consume an M operator
                        mx += dl, my += dl, ml -= dl;
                    }
                } else if (op == 'I') {
                    sy += SEQ.substr(cy, cl), sx += Array(cl+1).join("-");
                    cy += cl, my += cl, ++k; // consume an insertion I to ref
                } else if (op == 'S') {
                    cy += cl, my += cl, ++k; // consume a soft-clipping S
                } else throw Error("inconsistent MD"); // not working with N
            }
            if (ml != 0) throw Error("inconsistent MD");
        }
    }
    if (cx != mx || cy != my) throw Error("inconsistent MD");
    return [sx, sy];
}

var cigar = '16025H3M2I55M1D108M306H';
var MD = '7T15C24T9^C0G2T6T16T3A1A7G21T17G6C1T6T5A4';
var seq = 'GATCACAGGCCTATCACCCTATTAATCACTCACGGGAGCTCTCCATGCATCTGGTATTTTTTCGGGGGGGGATGCACGCGATAGCATCGCGGGCCGCTGGAACCGGAGCACCCTATGTCGCAGGATCTGTCTTTGATTCCTACCTCATGCCATTATTAATCGCGCCTA';
var s = md2pad(cigar, MD, seq);
console.log(s[0]);
console.log(s[1]);

great reference code with the caveat that MD is indeed really annoying to work with. oftentimes easier to just recreate the mismatches from the sequence itself

This program in C++ is able to perform the task https://github.com/mlafave/sam2pairwise

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


text:
How to interpret the SMRT Link base modification algorithm output?
I have got PacBio data from a Sequel machine at a coverage > 50X. I have run it in SMRT Link with their Base Modification Detection algorithm. I am now trying to figure out how to interpret the output and where to find the information about the modified bases. One of the outputs is a GFF3 file. The top of this file looks like this:
> ##gff-version 3

...
> I kinModCall  modified_base   1168    1168    31  -   .   coverage=48;context=TTAAACTAATTTCCAAAAATTTTGAACAGTTTTTAATACAA;IPDRatio=2.96
> I kinModCall  modified_base   1497    1497    32  -   .   coverage=50;context=CTTGAAAACACGAATATCTGTTTAGACAATAGTTTTGCCAA;IPDRatio=3.19
> I kinModCall  modified_base   1509    1509    33  -   .   coverage=48;context=CATGAAAAATTTCTTGAAAACACGAATATCTGTTTAGACAA;IPDRatio=2.84
> I kinModCall  modified_base   1529    1529    39  -   .   coverage=46;context=AAGAAAAGACCAAGAAAAACCATGAAAAATTTCTTGAAAAC;IPDRatio=3.84
> I kinModCall  modified_base   1904    1904    36  -   .   coverage=54;context=TAAAACCTGACAGAAAGAATTGGTCCCGTCATCTACTTTCT;IPDRatio=2.86
> I kinModCall  modified_base   2152    2152    31  -   .   coverage=59;context=TTAAAAATTAATATAAAATCAGGGAAATTTTTTTAAATTTT;IPDRatio=2.98
> I kinModCall  modified_base   2226    2226    36  -   .   coverage=63;context=TGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAAAA;IPDRatio=2.59
> I kinModCall  modified_base   2228    2228    34  -   .   coverage=52;context=CGTGGGCGGGGCGTCTCTGCCGGGAAAATTTTGCGTTTAAA;IPDRatio=2.98
> I kinModCall  modified_base   2237    2237    33  -   .   coverage=60;context=GCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAATTTT;IPDRatio=3.07
> I kinModCall  modified_base   2240    2240    44  -   .   coverage=63;context=GCGGCACGGTTTCGTGGGCGGGGCGTCTCTGCCGGGAAAAT;IPDRatio=2.89
> I kinModCall  modified_base   2311    2311    31  -   .   coverage=63;context=AAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAGAA;IPDRatio=4.35
> I kinModCall  modified_base   2313    2313    39  -   .   coverage=62;context=GAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATAG;IPDRatio=4.00
> I kinModCall  modified_base   2314    2314    34  -   .   coverage=63;context=AGAAAAATCGACGAAAATCGGTATAAAATCAAATAAAAATA;IPDRatio=2.63
> I kinModCall  modified_base   2333    2333    35  -   .   coverage=62;context=AAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGGT;IPDRatio=3.77
> I kinModCall  modified_base   2334    2334    35  -   .   coverage=63;context=CAAAAAAAAGAGAAAAAATGAGAAAAATCGACGAAAATCGG;IPDRatio=2.77

I expected something along the lines of 'chromosome I position 1497 methylcytosine', etc.
So it seems that the algorithm has detected modified bases at the following positions but it cannot call them, maybe because of a lack of coverage?
Update: the screenshot below shows there is no CSV file in the Data tab of the completed analysis.

It sounds like you're looking for the modifications.csv file, rather than the modifications.gff file. The .gff file is for viewing in SMRT View (or other viewers). You can see schema for both files at PacBio's Methylome Analysis Technical Note. Specifically, modified_base means the algorithm could not determine the exact base and modification.
I also found a video tutorial for performing the analysis in SMRTLink on YouTube
EDIT: Per PacBio's Methylome Analysis Technical Note documentation, the .csv "...is also produced when motif analysis is not active, such as with the RS_Modification_Detection analysis protocol." Not sure why it's missing, perhaps because the algorithm couldn't clearly identify any of the modifications. You'll probably need to contact tech support on this issue.
In either case, the .gff you posted includes the significant hits (p < 0.01), but the algorithm was unable to classify them, likely because of insufficient coverage. PacBio's Methylome Analysis Technical Note documentation gives specific coverage recommendations for each event type, but 100x is the lower bound. Lots of good information there.
Probably not what you wanted to hear, but I hope that gets you going.

I have run the Base Modification algorithm but I did not find the CSV file in the Data. Is that normal? I found all other files, including the gff file. I am trying again with the Base Modification and Motif Analysis algorithm. The analysis run was marked as successful.

@charlesdarwin, just want to make sure you saw my edit. Best.


text:
Aggregate sequencing/mapping/etc. metrics from cellranger across Illumina samples
I have a number of single-cell projects processed with cellranger from 10x Genomics. The pipeline produces a number of handy metrics that are summarized for each Illumina sample in web page. These include things like the estimated number of cells, the fraction of reads in cells, and mapping rates. What I would like is to aggregate this information from multiple Illumina samples into a single file. MultiQC is the obvious choice for such a task, but it appears that no one has written a module to process cellranger data for it. Are there any other prepackaged ways to perform this aggregation?
If there's no convenient tool for this I'll just write a MultiQC module, but I'd obviously like to avoid the extra effort unless it's really needed.

Single-cell RNA-seq analysis is still developing (not crystallized yet), so I think what you want is not there yet. If you think it may be handy for your analysis, many more may follow, so making a tool for this might be fruitful (if you make it available for others).

We’ve written a MultiQC module for Supernova, but not Long ranger or Cell Ranger (yet). I agree that one would be useful, and probably quite quick to write if there’s also a JSON summary file output for that tool too. We’ve been talking about this recently. Did you ever get anywhere with writing a new module?

any new about that? The cellranger is now version 5 and produce still data that would be great to compile in a single report.

Things move fast.  Cellranger now has a version 6.0

I wrote a simple script to do it. Not sure whether that's what you want.
#!/bin/bash
#CellrangerSummary.sh Sample1 Sample2
csvtool transpose $1/outs/*summary.csv | csvcut -c 1 > summary.csv
for sample in "$@"
    do
        csvtool transpose $sample/outs/*summary.csv > tempfilesforcsvtool.$sample
        csvjoin -c 1 summary.csv tempfilesforcsvtool.$sample > temp2
    mv temp2 summary.csv
rm tempfilesforcsvtool.$sample
    done

csvformat -T summary.csv > summary.tsv
rm summary.csv
```


text:
aligner for 1D^2 oxford nanopore data
I am wondering if anyone has tried aligning 1d^2 RNA-seq data. I am trying that on minimap2. It is giving me very low ~50% of mapped reads. it gives ~95% mapped reads for 1D.
I am using the following command for 1d and 1d^2:
minimap2 -ax splice -uf -k14 ref.fa direct-rna.fq > aln.sam

Update:
I also tried to play around the command using :
minimap2 -ax splice -k14 -a ref.fa direct-rna.fq > aln.sam
minimap2 -ax splice -ur -k14 ref.fa direct-rna.fq > aln.sam
minimap2 -ax splice -ub -k14 ref.fa direct-rna.fq > aln.sam

It did not work better.
FYI: I am trying to find answer for this from few days. I have posted this question on Biostars before. I did not get any response. Hence, I am trying to find the answer here.
Thank you in advance. I really appreciate any help.

Have you tried blatting/blasting a few of the non-aligning reads to check if they're okay?

No, I have not. But it will be good exercise. I will do it. Thanks.

You might consider trying to put the blast results into MEGAN to see if there is contamination in your sample. Also, could you please update us with what your sample type is (organism, tissue, origin), and if you have an amplicon, gDNA, cDNA library? Here is the link to MEGAN: https://ab.inf.uni-tuebingen.de/software/megan6

I have a similar problem (many unaligned reads). Did you find a solutioin? By any chance, did you compare minimap2 results with graphmap?

Please consider the answer @KVC_bioinfo and if you're happy with it for upvoting/accepting. Its good for site stats.

Since this is an old unanswered question where a lot of method development has since happened, going to put some answers here:

This rather old tool also mentioned in comments (graphmap) is supposedly good for this, though it may not be relevant anymore.
This paper uses minimap (not minimap2) alignments as a heuristic for read quality, suggesting that minimap2 should be good enough.

Note that depending on how you have processed (specifically basecalled) the data, 1d2 "second" reads may have low accuracy. That matches your ~50% unaligned proportion. You may need to do additional work during basecalling to fully use 1d2, as discussed in this paper.
Note: I am not an expert on nanopore and I put this answer together using google. Others are welcome to edit/editorialize.


text:
Can pair end reads with high heterozygosity be used to polish PacBio assembly?
I have used canu (correct) and smartdenovo to assembly PacBio reads. Next I am going to polish my assembly using illumina pair end reads by Pilon. However, I found there is high heterozygosity in my pair end reads when I used jellyfish to do k-mer analysis. There are two peaks in the k-mer frequencies distribution plot. The size of my genome is about 190M, and it is a diploid organism.

What I am doubting is that can I use this pair end data to polish my assembly?
If not, I have to do sequencing again.

In general polishing using short reads is depreciated because of mapping problems. If you have sufficient coverage, I would just stick to arrow.

Resequencing your sample is not going to make it any less heterozygous. I am assuming that you are working with a diploid organism? Please update your question and specify the genome size and ploidy of your organism (the type of organism will help us update our answers, too).
The answer, yes, you can use Pilon + short reads to correct the PacBio assembly even with a highly heterozygous sample. To do so you should use the --diploid option in Pilon to let it know that there will be heterozygous sites. The way that it corrects heterozygous SNPs and heterozygous indels is different, and it seems to perform best with SNPs and sometimes does not correctly identify heterozygous indels (source).
Pilon itself doesn't handle phasing, but just corrects with whatever reads are in the bam file that it is provided with. The end result of providing a bam with reads from both haplotypes is that the corrected sequence will be a mosaic of the dominant allele at the heterozygous site as was observed in the short read data. This may or may not select for some alleles based on sequencing coverage bias of Illumina sequencers.


text:
Can I stop my nanopore sequencing run if there are no more reads being produced?
I am sequencing a whole genome on MinION. I have used this flow cell for 6 hours for the phage lambda control. That gave me around 300,000 reads. I started a 48h sequencing run on the same flow cell after washing it. The sequencing has been going on for 1 day and 21 hours. It is not producing any reads at the moment (Strand = 0). The number of single pores is also 0. Can I "Stop acquisition" without negatively affecting the results?
The image below is a screenshot taken from MinKNOW, the software used to control the MinION.

Yes, you can click "stop acquisition" and your run won't be negatively affected. All of the reads are saved as they are generated. I am not sure how this will impact live basecalling though if that is something that you do.

Live basecalling for the uncalled reads will continue until all have been called, or the application is quit.

Yes, you can stop the run and it should not affect the base calling. We normally do the runs and then do the base calling after the run has completed.


text:
Does the DNA or RNA of some ethnic groups map better than others' to the human reference genome sequence?
I believe that the human genome reference took DNA samples from different people and that some natural variation is included in extra contigs.
However, the human reference genome comes from a limited number of people.
Can the limited number of donors and the linear representation of the genome (as opposed to a graph representation that could potentially include variation) bias an analysis in favour of one ethnic group over another? Are there studies that have compared the alignment rates/mappabilities of reads of people coming from different geographic regions? Do some ethnic groups' reads align better to the reference?

For simple variants like SNPs it would not really be a problem to use the current genome assembly for other ethnic groups. But for more complex variants this could be indeed problematic, however not only for ethnic groups but also for individuals within the same population. Think of very complex regions such as HLA or KIR.
In studies where they compare different ethnic groups (e.g., 1000 G) they rather not discuss such matters, but act like the analysis is 100% perfect.
There is a nice comment in Genome Biology about this issue.


text:
Impute phenotype under some constraints
Errors happen frequently in the lab and I got one sample that we mixed one important information (the region of the sample) between two samples. Now we don't know which sample is from which region, but we know all the other variables (Date of the sample, patient, and other variables).
So far, my solution has been to make that values as NA to be able to analyze it, but I would like to know if there is any solution. If I can convey to test or other analysis downstream that sample A and sample B are from Location1 or Location2 and if A is from Location1, B is from Location2 and otherwise round.
There are several packages for imputing quantitative data, but I couldn't find any to store this kind of information or to impute qualitative data. I can recall a recent paper with single-cell sequencing that estimated phenotype and some other variables but I can't find it neither I think it is possible with RNAseq.
Although maybe I could build a classifier to actually do that (but that would be too long and I'm not sure if it is worth it)
Notes

My samples are from more than two locations.
My samples consist of 16S microbiome DNA-seq and RNA-seq

Question 1: Is there a tool for imputing phenotypes while also using genetic data?
Yes, try Andy Dahl's Phenix from the Marchini lab.
Disclosure: I did my PhD in the Marchini lab, although I was not directly involved in this research.
Question 2: Is there a way to impute categorical data using methods intended for continuous data?
If the variable of interest has only two categories, then you should be able to simply recode the categories as arbitrary (differing) numbers and get sensible results. But do make sure you understand the tool you are using to see if this will work.
For more than two categories, if your imputation approach allows the simultaneous imputation of multiple variables, you could try recoding your categories as contrasts or one-hot.  After recoding, you could impute the variables resulting from recoding together.


text:
Can I pilon-polish long reads with Illumina short reads to improve structural variant detection?
I have pacbio Sequel data at 50 x coverage for a strain of animal. I would like to find structural variants compared to the reference genome sequence. At the moment, I align my reads to the reference genome using NGM-LR and call the variants using Sniffles. The resulting VCF file contains variants that are well supported in a visual inspection in IGV as well as some more questionable variants (especially larger events). I have used the reads corrected by canu for alignment. I think it splits my long reads into very long reads and not-so-long reads and uses the shorter reads to correct the longer reads. However, in IGV, my reads still have a considerable number of mistakes compared to the reference. I was wondering if it is OK to use pilon to polish my long reads with some Illumina data from the same sample? The main purpose of pilon is to correct assemblies but I was wondering if it can be used to correct subreads from PacBio sequencing?

I have no links to support it, so I will give just a comment. People say that using short reads to correct have always troubles with misplacing them due to repetitive content. I suspect that structural variants could be frequently duplicated and therefore problematic for error correction with short reads and I expect that with 50x you will get better error correction using just the long reads.

I'd recommend using Canu for correction rather than Pilon, because it has a component that is specifically designed for read-level correction. The newest version of Canu (v1.7) will track all reads through the correction and assembly process:

Ensure that every raw read is either corrected or used as evidence for correcting some other raw read. This serves to rescue short plasmids in high coverage datasets, and it is no longer necessary to set corOutCoverage to achieve the same result.

More details here.


text:
Extracting the CIGAR string from a BAM via Python?
Is there a standard method in Python to extract a CIGAR string from the BAM?
There are great libraries which parse the CIGAR, e.g. https://pypi.python.org/pypi/cigar/0.1
>>> c = Cigar('10M20S10M')
>>> c.mask_left(10).cigar
'10S20S10M'
>>> c.mask_left(9).cigar
'9S1M20S10M'
>>> Cigar('10S').mask_left(10).cigar
'10S'
>>> Cigar('10H').mask_left(10).cigar
'10H'
>>> Cigar('10H').mask_left(11).cigar
'10H'
>>> Cigar('10H').mask_left(9).cigar
'10H'

It looks like pysam already parses the cigar string, if I'm not mistaken:
import pysam

bam = 'myfile.bam'
bamfile = pysam.AlignmentFile(bam, 'rb')

for read in bamfile:
    if not read.is_unmapped:
        cigar = read.cigar
        print(cigar)

This outputs lists of tuples of the parse output:
[(5, 61), (0, 30), (5, 198)]
[(4, 11), (0, 30), (4, 248)]
[(4, 11), (0, 30), (4, 248)]
....

How could I simply output the CIGAR from read? I would prefer to create a column in a pandas DataFrame

If you really do just want the cigar string then it's read.cigarstring. However, I'm not sure what you're trying to gain with the cigar package from Brent. Unless you want to get the string with the masking changed then the tuple you get from pysam is the same as what you get from cigar (with the exception of the numeric operations instead of character operations).

Oh, the API on the website had me confused with cigar string. Perhaps I could use the parsed cigar from pysam, but I'm confused: How do I parse the list of tuples in order count the total number of insertions, deletions,  matches (and possible other information) from this?  I apologise if this is documented, but the pysam manual is confusing imho.

Should I read the above, e.g. the first line [(5, 61), (0, 30), (5, 198)] as 61 hard-clippings, 30 matches, and 198 hard clippings?

Yes, that is the correct interpretation of that list of tuples, which is documented in the API

I see; this question is largely my misunderstanding of the pysam website. I'm happy to delete this question, unless it's helpful for the community.

# pip install pyranges
# pip install bamread

# wget http://hgdownload-test.cse.ucsc.edu/goldenPath/hg19/encodeDCC/wgEncodeRikenCage/releaseLatest/wgEncodeRikenCageGm12878CellPapAlnRep1.bam # (half a gig)
# samtools index wgEncodeRikenCageGm12878CellPapAlnRep1.bam

f = "wgEncodeRikenCageGm12878CellPapAlnRep1.bam"

import pyranges as pr
df = pr.read_bam(f, sparse=False, as_df=True)

df[["Chromosome", "Start", "End", "Strand", "Cigar"]].drop_duplicates("Cigar")
#          Chromosome  Start    End Strand             Cigar
# 0              chr1  10579  10606      -               27M
# 36             chr1  29334  29360      +           8M1I18M
# 37             chr1  29334  29360      +           7M1I19M
# 38             chr1  29334  29360      +           5M1I21M
# 44             chr1  29339  29365      +           4M1I22M
# ...             ...    ...    ...    ...               ...
# 19655298       chrM  14550  14575      +       16M5I3M3D3M
# 19660213       chrM  14747  14775      -       15M5D5M4I3M
# 19663987       chrM  15374  15399      -   2M1I2M2D4M3I15M
# 19676864       chrM  16234  16271      +  15M1I4M12D5M1I1M
# 19677003       chrM  16264  16291      +   14M1D3M2D1M3I6M


text:
seqret error when trying to transform GFF3 and FASTA files into EMBL format
I am using the tool seqret from emboss to transform an annotation file in GFF3 format and a FASTA file into an EMBL file because Wormbase does not supply an EMBL file with annotation and sequence together. I am running it on Mac OS X. I got the gff3 and fasta files from Wormbase. I installed emboss using conda from the bioconda channel. I use the following command:
seqret -sequence c_elegans.PRJNA13758.WS263.genomic.fa -feature -fformat gff -fopenfile c_elegans.PRJNA13758.WS263.annotations.gff3 -osformat embl -auto -outseq c_elegans.PRJNA13758.WS263.embl
It fails with the following error:
/Users/me/anaconda/bin/seqret`: line 7: 54774 Killed: 9               _seqret "$@"

/Users/me/anaconda/bin/seqret is a short file that indeed contains _seqret "$@" on line 7:
#!/bin/sh
BIN_DIR=$(dirname $(which seqret))
export EMBOSS_ACDROOT=$BIN_DIR/../share/EMBOSS/acd/
export EMBOSS_DATA=$BIN_DIR/../share/EMBOSS/data/
export PLPLOT_LIB=$BIN_DIR/../share/EMBOSS/
_seqret "$@"

I tried to erase it but that did not help.
I also tried EMBLmyGFF3 in a Python 2 environment:
EMBLmyGFF3 c_elegans.PRJNA13758.WS263.annotations.gff3 c_elegans.PRJNA13758.WS263.genomic.fa --topology linear --molecule_type 'genomic DNA' --transl_table 1 --species 'Caenorhabditis elegans' --locus_tag CELE --project_id PRJNA13758 -o c_elegans.PRJNA13758.WS263.annotations.embl
but I just had this message on display for hours:
#############################################################################
# NBIS 2016 - Sweden                                                        #
# Authors: Martin Norling, Niclas Jareborg, Jacques Dainat                  #
# Please visit https://github.com/NBISweden/EMBLmyGFF3 for more information #
#############################################################################

Update 11th March:
I started EMBLmyGFF3 yesterday at 5 pm and this morning it started outputting messages such as:
08:22:08 ERROR feature: >>inverted_repeat<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>assembly_component<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'genbank' - skipped
08:22:08 ERROR feature: >>tandem_duplication<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'public_name' - skipped
08:22:08 WARNING feature: Unknown qualifier 'polymorphism' - skipped
08:22:08 WARNING feature: Unknown qualifier 'other_name' - skipped
08:22:08 WARNING feature: Unknown qualifier 'consequence' - skipped
08:22:08 WARNING feature: Unknown qualifier 'variation' - skipped
08:22:08 ERROR feature: >>biological_region<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 WARNING feature: Unknown qualifier 'balancer_type' - skipped
08:22:08 WARNING feature: Unknown qualifier 'balancer' - skipped
08:22:08 ERROR feature: >>protein_match<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>transcript_region<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:22:08 ERROR feature: >>transcribed_fragment<< is not a valid EMBL feature type. You can ignore this message if you don't need it.

...
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:45:36 ERROR feature: >>piRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:46:00 ERROR feature: >>pseudogenic_tRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:47:59 ERROR feature: >>antisense_RNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:47:59 ERROR feature: >>lincRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:48:03 ERROR feature: >>pre_miRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:48:03 ERROR feature: >>miRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
08:50:37 ERROR feature: >>snRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
09:04:00 ERROR feature: >>miRNA_primary_transcript<< is not a valid EMBL feature type. You can ignore this message if you don't need it.

For some of these features, it seems OK to not report them. But for others, such as piRNA, I need to find a solution. C. elegans has got 1000s of piRNAs, it would be nice to have them in the EMBL file. I need to find a list of valid EMBL feature types.
Update + ~ 24h computation:
...
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
12:21:19 ERROR feature: >>scRNA<< is not a valid EMBL feature type. You can ignore this message if you don't need it.
Otherwise tell me which EMBL feature it corresponds to by adding the information within the json mapping file.
Conversion done

My EMBL file is only 1.19 GB compared to the GFF3 3.34 GB, so I probably lost a lot of annotation that is not supported by EMBL.

Please always mention your operating system and hardware. I am guessing that you are running this on a mac laptop, but it would be much easier if I didn't need to guess.

That error is telling you that the script you were running was killed by signal 9, also known as SIGKILL, while executing the command _seqret "$@". This is the way the operating system will kill processes that are asking for more memory than it can provide.
Note that removing _secret "$@" is not the solution since that's what's actually doing the job. $@ is a bash array that contains all the arguments passed to a bash script. The /Users/me/anaconda/bin/seqret file is just a simple script that sets up a few variables and then launches _seqret, passing it the arguments you gave to /Users/me/anaconda/bin/seqret as $@. So deleting that lines just makes the script useless.
Based on your question, I am guessing you're running macOS on a laptop, so it isn't unreasonable that something that needs to read an entire genome and what looks like a full set of annotations will run out of memory on a laptop. Of course, this is a pretty small genome, but the annotations are still 396MB compressed, so depending on your machine's specs and what kind of internal processing the command is doing, it isn't too surprising that it fails.
So, either run this on a more powerful machine, or try splitting the annotations file into smaller files and running it sequentially on each of those.

I'm one of the EMBLmyGFF3 developers. If the compressed annotation is 396MB as mentioned by @terdon it is not a small annotation... In the worse case we have tested we never waited more than few minutes. We haven't done any particular speed test. It would be valuable if you could tell us how long it has taken in your case (If you didn't kill the process before the end...).
At least we could try to implement some nice features, like to see that the tool is working , or/and see which percentage of the file has been already processed.
If it is really slow we could think to try to improve its speed.

Does EMBLmyGFF3 accept .gz compressed files? Is it better? I have started the transformation at 5 pm today. The size of the uncompressed annotation file is 3.1 GB. The genome file is about 100 MB uncompressed. I'll let you know if/when it finishes. It would definitely be useful to have a progress bar and a sign that the tool is working properly. Thanks a lot!

Still going and using 50 GB of memory (on a 16 GB physical RAM machine) but less than 100% of a CPU. So, it's using the hard disk to store some stuff.

The process has taken about 24 hours in total on 16 GB physical RAM and 8 cores.

Ok it is obviously really slow... Using .gz compressed file will not change anything here. We will investigate how to improve the speedness.

I get a lot of errors of the type >>antisense_RNA<< is not a valid EMBL feature type. I found http://www.insdc.org/files/feature_table.html#2.3 but I could not locate a table with all the valid EMBL feature types. Do you know where I could find this? I get the error with miRNA as well, would ncRNA be the correct feature type in the EMBL format?

You will find the complete list of features and qualifiers accepted by EBI here: https://www.ebi.ac.uk/ena/WebFeat/


text:
Intuitive explanation of GSVA analysis
I'm trying to understand the way the GSVA analysis is working behind the scenes.And I was wondering if there is any way to understand it more intuitively the whole process.
So at first according to paper it starts by evaluating whether a gene i is highly or lowly expressed in sample j in the context of the sample population distribution. They use these kernel estimations of the cumulative density functions to transform the initial values so not to be affected by the problematic intensities.
After this "transformation" and a following normalization, GSVA calculates the enrichment scores using the Kolmogorov-Smirnov (KS) like random walk statistic.
As I know, Kolmogorov-Smirnov checks for differences in distributions. Which distributions does it check? Gene-set's against all the others genes? And what is the role of the random walk?
So is there any intuitive way to understand this kind of Kolmogorov-Smirnov (KS) like random walk statistic? How does it actually work? Which one is the null and which the alternative hypothesis in that case?

The most intuitive explanation is  also explained in the background section:

Conceptually, this methodology can be understood as a change in coordinate systems for gene expression data, from genes to gene sets.

Which I think it is explained in Figure 1 of the paper, which is also one of the most informative I found about the methods of an algorithm1.

I'll will use it to explain what's going on.

Fit a distribution function per gene

For each gene expression profile ... a non-parametric kernel estimation ... function is performed.

I can't give more detail because I haven't understood it fully

Then they are normalized to make the ranks symmetric around zero
It is an inline formula between equations 2 and 3.
The Kolmogorov-Smirnov like random walk statistic is applied to those kernel estimations normalized distributions.
Which two distributions does it compare? The genes in a gene set and all the others (but I'm not totally sure)
The role of random walk is explained a bit later from your quote:

[the random walk] produces a distribution over the genes to assess if the genes in the gene set are more likely to be found at either tail of the rank distribution

Which is later used to calculate the enrichment score by looking the maximum deviation from zero, or by summing the largest and the lowest deviations.

To address some other questions:

GSVA calculates sample-wise gene set enrichment scores as a function of genes inside and outside the gene set, analogously to a competitive gene set test

1: Try finding which fittings and approximations are done inside limma


text:
Evaluate clusters of individuals by using their sequence data
For a dataset of several hundred individuals, I applied a hierarchical clustering to generate clusters based on a functional trait that sets them apart.
My task is now to evaluate if these clusters can be supported by the nucleotide  sequence data of the corresponding gene, i.e. is there more genetic similarity within each cluster than between the clusters.
For a first approach, I created a multiple sequence alignment for every cluster and calculated the the % identical sites value. There are more identical sites within each cluster than in an alignment of all sequences.
Do you know any good tools (preferably python) that can perform a more sophisticated evaluation of existing clusters based on sequence data?
Are there other scoring methods with which I can evaluate the sequence similarity within the clusters (again preferably python)?

You can calculate allele frequencies for each cluster you have to further verify if they belong to similar population, however if size of your dataset is rather small this may not work for you.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3813878/, here is an article about differentiating populations based on sequencing data.


text:
Convert bam file to highly compressible bam
I have a large collection of bam files and I want to post-process each of them into another bam where I can make queries about:

the reads position and pair-endness,
insert sizes,
MAPQ and other flags, etc. of the reads

but where I don't need to make use of:

the sequence string of the reads,
the Q-scores qualities or,
the read ids, apart from keeping pair-endness in the read id,

So I would like to know what's a good tool/parameters to process each bam and turn them into a small size bam that still retains the information needed above. CRAM format is not an option at this point as several of the tools I use are not CRAM-compatible yet.
Any ideas?

Have you considered third-party BAM compression solutions?

You can just set the fields you don't need to *:
samtools view -h foo.bam | awk 'BEGIN{FS="\t"; OFS="\t"}{if($1~/^@/) {print $0} else {print "*", $2, $3, $4, $5, $6, $7, $8, $9, "*", "*"}}' | samtools view -bo smaller.bam -

This will set each read's name to *, but you can still see where its mate maps with the PNEXT and RNEXT fields.
The resulting BAM file looks like:
*   83  chr1    3003308 37  51M =   3003157 -202    *   *
*   99  chr1    3003364 33  51M =   3003547 234 *   *
*   99  chr1    3003410 37  50M =   3003636 269 *   *
*   99  chr1    3003438 255 50M =   3003760 373 *   *

Why *, specifically? Is that part of the bam specs or will any character do?

* is defined in the BAM specification to mean "unknown" or "not available". You can get similar results if you convert to CRAM and don't decode the read name or seq/qual fields (i.e., you'll get a * in those fields).

I wrote  http://lindenb.github.io/jvarkit/Biostar173114.html for biostars "What is the best way to make a bam file smaller by removing unwanted information?"
Usage: biostar173114 [options] Files

    -keepAtt, --keepAttributes
      keep Attributes
      Default: false
    -keepCigar, --keepCigar
      keep cigar : don't remove hard clip
      Default: false
    -keepName, --keepName, --name
      keep Read Name, do not try to create a shorter name
      Default: false
    -keepQuals, --keepQualities
      keep base qualities
      Default: false
    -keepRG, --keepReadGroup
      if attributes are removed, keep the RG
      Default: false
    -keepSeq, --keepSequence
      keep read sequence
      Default: false
    -mate, --mate
      keep Mate/Paired Information

e.g:
$ java -jar dist/biostar173114.jar src/test/resources/S1.bam --keepSequence | head -n 20
@HD VN:1.5  SO:coordinate
(...)
@SQ SN:RF11 LN:666
@RG ID:S1   SM:S1   LB:L1   CN:Nantes
R0  128 RF01    1   60  70M *   0   0   GGCTATTAAAGCTATACAATGGGGCCGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACT  *
R1  64  RF01    8   60  70M *   0   0   AAAGCTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATG  *
R2  64  RF01    11  60  70M *   0   0   GCTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATCCGC  *
R3  64  RF01    12  60  70M *   0   0   CTATACAATGGGGAAGTATAATCTAATCTTGTCAGAATATTTATCATTTATATATAACTCACAATCCGCA  *
R4  128 RF01    27  60  70M *   0   0   GTATCATCTAATCTTGTCATAATATTTATCATATATATATAACTCACAATCCGCAGTTCAAATTCCAATA  *
R5  64  RF01    44  60  70M *   0   0   CAGAATATTTATCATTTATATATAACTCAGAATCCGCAGTTCAAATTCCAATATACTATTCTTCCAATAG  *
R6  128 RF01    44  60  70M *   0   0   CAGAATATTTATCATTTATATATAACTCACAATCCGCAGTTCAAATTCCAATATACTATTCTTCCAATAG  *

Thanks @Pierre, what happens if it's run without the --mate flag?

the mate is included (default is : set read a single read)


text:
Salmon output: transcripts quantified with zero reads support
I quantified some samples using Salmon. According to the documentation of output format, the last column of the 'Quantification File' represents the number of reads that supported the given transcript, namely:

NumReads — This is salmon’s estimate of the number of reads mapping to each transcript that was quantified. It is an “estimate” insofar as it is the expected number of reads that have originated from each transcript given the structure of the uniquely mapping and multi-mapping reads and the relative abundance estimates for each transcript.

In some case, I had in the output some transcripts that had NumReads=0.000000. Does this mean that I should discard them from the list of transcripts obtained in that sample? How comes that those transcripts end up in the output file in first place, since they don't have any reads support and the quantification was done using the reads themselves?

All transcripts should always be returned in the output, whether there's evidence supporting their expression or not. This is for the sake of simple convenience, since typically one would be merging these values across samples into a matrix. Missing values will then tend to become either NA or cause an error, when what would be appropriate would be to have values of 0. This also facilitates calculations such as "what's the median number of reads per transcript?" or "what percent of transcripts have an expected reads count of at least X".

Ok thanks. Only the transcripts with reads support should be then consider as 'expressed' right? And what about read support value between 0 and 1? Shall I consider them to be present or not?

Normally one would use something like zFPKM to come up with some "expressed" threshold. Note that values above 1 don't necessarily mean that there is even a single read uniquely assignable to a given transcript.


text:
While opening a Bam file with the SeqAn library I get 'seqan::FileOpenError'
When I try to open a simple bam file with the SeqAn library using the following code:
CharString bamFileInName = getAbsolutePath("/home/omair/TFG/Files/9827_2#49.bam")
BamFileIn bamFileIn(toCString(bamFileInName));

I get the following error:
terminate called after throwing an instance of 'seqan::FileOpenError'
  what():  Could not open file /home/omair/Programs/seqan//home/omair/TFG/Files/9827_2#49.bam: iostream error

Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

The file exists and the path is correct.

cross posted: https://www.biostars.org/p/302804/

Do you really have a file called /home/omair/TFG/Files/9827_2#49.bam? Please [edit] your question and show the exact code you are using. What you have posted has to be a copy/paste error, it makes no sense.

Duplicate of https://bioinformatics.stackexchange.com/questions/3804/cannot-open-bam-files-with-the-seqan-library

Remove the getAbsolutePath call — your path is already absolute. getAbsolutePath garbles it, as you can see in the error message.

By removing the getAbsolutePath call I get the following error:   terminate called after throwing an instance of 'seqan::UnknownExtensionError'
what():  Unknown file extension of /home/omair/TFG/Files/9827_2#49.bam: iostream error
Process finished with exit code 134 (interrupted by signal 6: SIGABRT)

@OmairIqbal Ensure that seqan was compiled with zlib support. See here for details.


text:
simple data subsetting on r
I have a big df like following (pseudo) df and I want to retrieve values in different col. according to the rows in the first col.
A B C D
a 1 1 1
a 2 3 1
a 3 3 3
a 4 5 6
b 5 6 5
b 7 6 5
c 7 6 5
c 5 4 3
c 5 5 4

I want to retrieve the values in col. B, C, and D according to the rows in col A and perform further analysis of them. For example, multiply the values in col B and C for all the rows of "a" in col A. Can someone suggest a commend for such a naive question. Thanks.

I'm downvoting this because it doesn't have anything to do with bioinformatics, but you can use df[df$A=='a',], or dplyr::group_by(A) if you want to repeat the same analysis for each value of A

That is not true, the original df contains genetic info. Thanks

I'm voting to close this question as off-topic because without more context, this looks like a programming question better suited for Stack Overflow, and the "this question belongs on another site in the SE network" radio button does not offer Stack Overflow as a selection.

Here's an approach that uses native R, using the vectorized (fast) family of apply() functions.
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> apply(df2, 1, prod)

Use apply/lapply/mapply/etc. where you can in place of alternatives, where you don't know their run characteristics.
Here's a more concrete example:
> df
  A B C D
1 a 1 1 1
2 a 2 3 1
3 a 3 3 3
4 a 4 5 6
5 b 5 6 5
...
> df2 <- df[df$A == 'a', which(names(df) %in% c("B","C"))]
> df2
  B C
1 1 1
2 2 3
3 3 3
4 4 5
> apply(df2, 1, prod)
 1  2  3  4
 1  6  9 20

Edit: If you need to do this on each unique value of df$A, you can nest apply calls to make this scale:
> perA.list <- apply(data.frame(unique(df$A)), 1, function(x) apply(df[df$A == x, which(names(df) %in% c("B","C"))], 1, prod))
> perA.list
[[1]]
 1  2  3  4
 1  6  9 20

[[2]]
 5  6
30 35

[[3]]
 7  8  9
35 20 25

Each index in this list correspondents to the index of the value from unique(df$A), i.e. a, b, and c.
If you don't want to do this for all unique values in df$A, but for some specific element in df$A, then you would just need to write a function. 95% of the code is in this answer to write such a function, so I'll leave that as a homework exercise for you to implement.

Thanks, I can't able to mention "a"every time. First I also want to have same results for "b" and "c" as for "a". Second, as I mentioned the data is big with 70,000 entries. In this case please suggest something.

See edited answer.

Hi, ddply (https://www.rdocumentation.org/packages/plyr/versions/1.8.4/topics/ddply) was used to obtain results, exactly what I was looking. Thanks


text:
Command not found error on HISAT2
I have been attempting to map reads to a reference genome using HISAT2 using the Pertea, et al 2016 Nature Protocols paper. I am able to successfully generate an index (using the hisat2-build command), but when it comes to the actual mapping, I run into problems. The code I'm using (from the paper):
$hisat2 -p 8 --dta -x daphnia_tran -1 C2*R1*.fastq -2 C2*R2*.fastq -S C2paired.sam

And I get the error message:
./HISAT2_Basic_14.sh: line 26: -p: command not found

When I remove the -p argument, the error message is similar, but says "--dta: command not found"
Any ideas on what the issue is?
Thank you!

tldr: Remove the $ from the command.
I imagine you're literally typing $hisat2, where you mean to instead type hisat2. The $ is just meant to show the end of the command prompt. If you actually type $hisat2, you're not executing the hisat2 program, but whatever the $hisat2 variable is set to. Since it's likely not set to anything, it's just ignored. For example:
ryan@maximus:~$ $hisat2 ls
bin              convertChromosomesChain.py  ...

Since there's no $hisat2 variable, I'm actually just running ls. Likewise, your command is just reducing to -p 8 --dta -x daphnia_tran -1 C2*R1*.fastq -2 C2*R2*.fastq -S C2paired.sam, which is the source of the error.

Can't believe I missed that, thanks very much! The script seems to be running perfectly now.


text:
Two variants associated with same chromosomal position in genotype data for one individual
I'm trying to play around with my family's raw genetic data from Ancestry.com so that I can do some genetic analyses, and I noticed in the AncestryDNA raw data files that some regions of the genome have more than one variant associated with them and I would like to know why this is. The end goal is get this into a regularly formatted VCF, but I feel that these double variants will probably mess things up. So for example one might see something similar to this in the tab-delimited file:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

(Header line: rsid    chromosome    position    allele1    allele2)
I understand that when rsid numbers were created, sometimes multiple rsid numbers were accidentally assigned to the same region of the genome, but  I do not understand how two separate variants can be assigned to one single region in the genome. I could be wrong but as I understand it, I should only see one line here (either a "T    T" or an "I    I") because genotypes are reported in comparison to the human reference genome and so each position should uniquely refer to a single nucleotide in the said human reference genome.
But in this example case one line says there is a "T    T" genotype at this position, and the other line says there is an insertion(see "How do I read my raw DNA data?" here).  One possible explanation I thought of is if this would just be saying that a T was inserted, but if this is the case, then as I understand it the reference allele present in the individual should still reported, just in conjunction with the allele at the index. (So assuming an extra "A" was inserted after the T, which could be the reference or alternative allele, something like this should be recorded: rs6010164  22  50515270    TA  TA).
Can anyone help me understand why there are two separate variants associated with chromosomal position? Thanks!

Are you sure that if an A was inserted, you would get GA GA? I have seen G GA for an A insertion. Looking back at your header, it does not seem to respect the VCF format. The order of the first three columns is scrambled. It should be chromosome, position, id. Then, there should be the reference sequence at that position, followed by the alternative sequence. Then, in the sample column, you should get the genotype in function of the reference and alternative alleles.

Hi @charlesdarwin yes you are right its not a regular vcf, I should probably just refer to it as a tab-delimited file. No I am not sure that I would get that, I am not entirely sure what format indels are supposed to be in, but I do feel like its super weird for them to have more than one line associated with a single position, and eventually I need to change this into a VCF file and I feel like that would probably mess things up.

Did they perform exome sequencing or use a SNP chip? This is likely to occur for the latter, since they're then reporting every probe.

@DevonRyan I believe they used a SNP chip.

The VCF format does allow multiple variants at the same chromosome & position - in fact the official dbSNP VCFs have exactly that.

First of all, take these data with a very large pinch of salt. This sort of targeted analysis is not designed to produce high quality genetic data but, to give an idea of a sample's ancestry. Given a file like the one you describe with no quality information, you have no way of knowing how reliable any of the variants called actually is. From personal experience with another, similar company, I expect that they aren't reliable at all.
I have tested four individuals using a commercial tool like the one you describe and while I don't doubt the validity of their main results (those displayed on their website), their raw data showed combinations of variants that were extremely unlikely to be real given that the individuals were still breathing.
This isn't to say that such companies don't do their job well, only that their job is limited to what they actually display. When you look into their raw data, you're on your own. Presumably, they apply various filters before displaying the main results, but these are most likely not applied to the raw data. So if you want to do your own analyses, this is unlikely to be a reliable dataset. Also, I would strongly urge you to seek the help of an expert when interpreting any data you do have.

That said, the rsIDs are simply identifiers of a specific genetic variation, not of a location. Any genomic position can have many (essentially infinite) possible variations. For example, a C residue at position 100 of chromosome 1 could become a G, a T or an A. Each of these will have its own, separate rsID. Now, if you also consider insertions and deletions, you can see why you can have essentially limitless possible variations at that specific position. So no, having more than one rsID describing the same position is not strange at all.
The way these tests work is they look at a list of possible variants, each of which has its own rsID, and they report what they find for that rsID in the sample. Then, they will report what they find at the rsID's position for each allele. This is explained in the FAQ article you linked to:

The possible observations are A for adenine, C for cytosine, G for guanine, T for thymine, I for insertion’ and ‘D for deletion, or 0 for missing data. Column one provides the identifier (including the #rsID where possible). Columns two and three contain the chromosome and basepair position of the variant using human reference 37.1 coordinates. Columns four and five contain the two alleles observed at this variant (genotype)

Let's have a look at the two rsIds you mention.

rs6010164: This is a single nucleotide variant (SNV) which describes a C residue where the reference genome has a T.
rs267607236: this is a deletion at position 50076841, so a very different variant to the above. In this one, the sequence TGAG is deleted.

Your results are:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

So, your sample doesn't have rs6010164 since it has the same sequence at that position as the reference. However, it does have the insertion as indicated by the two Is.
The Is and Ds of the insertions and deletions (indels) are going to be the easy bit here. The only way you can parse this sort of file and convert it to VCF is to get the list of target SNV rsIDs, check what they describe as the reference and as the variant and then, for every one, check which of the two is shown in your file.
If you're not an expert in this field, you will probably be better off hiring someone to do this for you. There are actually companies that offer this sort of service and can take these tab-delimited files as input. I even work for one!

It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here.
Without info about the reference allele it's pretty confusing (particularly because they appear to be reporting the opposite strand from the database (G vs. C/T), but I would interpret this to mean that both alleles match the reference for rs267607236 (they have the insertion), but do not match the reference for rs6010164 (G-C pair rather than T-A). This is also the most common genotype at this position, being found in almost 80% of people in most populations.

Thank you @heathobrien! If I want to convert my data to a VCF don't I need to one use only one of the rsid's since they are at the same position? Is there any standardized way to pick which one to use? And sorry, my bad, they did report from the forward strand.  I had changed the nucleotides since its actually not my own DNA data I'm posting here, but I should have checked to make sure I was changing it to a genotype that made sense.

You need to understand a few key points for your example given:
(a) Ancestry (and 23andMe) will use the letters I and D to indicate it is an InDel/DELINS. You have to read the rsID definition to determine if the variant is an Insert or Delete from the reference.  And therefore which letter is used to mean positive / variant and which for negative / reference.  They simply use the other letter to indicate negative / reference.
rs267607236 is defined as a deletion. So the I's mean the reference was detected in this case. Therefore the Ancestry file indicates no deletion of the 4 base-pairs at coordinates 50515271-4.
(b) Ancestry and dbSNP VCF annotation files uses a "left-aligned" coordinate for ALL InDel's and generally designate the coordinate BEFORE the variant.
This makes sense for Inserts. It gives the coordinate of the last base-pair in the reference just before the insert starts. That coordinate plus 1 is the next base-pair after the insert in the reference.
This same method is used for deletions.  But because the deleted values have coordinates in the reference, this makes deletions appear "one off". You could use the starting coordinate of the deletion when defining the variant (e.g. [TGAG/-]). Otherwise, you need the coordinate before if defining the deletion as [CTGAG/C].
rs267607236 is listed in the OP's information as coordinate 50515270. Because it is a deletion, the deleted sequence actually starts with the base-pair at coordinate 50515270+1 or 50515271.
(c) As it ends up, for rs267607236 in question here, the dbSNP VCF and the dbSNP online website entry are different (both latest release 156 in this case). This could happen but still be consistent because there is a T flanking the 4 base-pair deletion.  The starting reference sequence ...C TGAGT A... becomes CTA after the deletion.  But which T remains is not really known. The VCF file defines the deletion by retaining the trailing T. The online website retains the leading T.
And in fact, it appears the website entry starting coordinate is incorrect.  Given the way they defined this variant within the flanking sequence, the online website should use the coordinate 22:50515270 (like in the Ancestry file) and not 22:50515271 as it does now. This may have existed with that coordinate in an earlier release used by Illumina / Ancestry when defining the pipeline.  But why the VCF and online website are different now is confusing.
There are numerous other nuances, conventions and actual errors that have to be understand to properly convert a microarray raw data file to a VCF.
Here are all the entries that we referenced / consulted brought together for your perusal:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

GRCh37 Reference sequence around rs267607236:
TTAC TGAGT AAGA  (70-74)

dbSNP Release 156 GRCh37 VCF:
NC_000022.10    50515270        rs6010164       T       C
NC_000022.10    50515269        rs267607236     CTGAG   C       # Note: leading T deleted

dbSNP web page definition:
Alleles:T>C             for rs6010164)
TTAC [T/C] GAGTAAGAT
22:50515270 (GRCh37)
NC_000022.10:50515270:T:C

Alleles:TGAGT>T         for rs267607236
TTACT [GAGT/-] AAGAT                                            # Note: trailing T deleted
22:50515271 (GRCh37)                # for consistency with flanking sequence, this should be ..270
NC_000022.10:50515271:TGAGT:T

It is not clear to me what Ancestry does when the same rsID has an insert and a deletion at the same position. For example, the line from dbSNP VCF for rs768019142 is:
NC_000001.11    10055   rs768019142     TA      T,TAA
which says there are two variants in the rsID. (This is common, especially for SNPs with multiple variants possible).  A deletion at 10056 OR a single letter insert after 10056. To represent this, instead of a single letter for the reference, they have given two base-pairs as the reference sequence. They probably should have defined two rsIDs at two adjacent locations. UCSC's versions of the VCF files are filtered and cleanup to simplify items like this.

For different variation type at the same location, each variation type is being assigned its own rsID by the dbSNP maintainers.
In your example, one chromosome in the pair has a single-nucleotide variation (T allele), while the other chromosome has an insertion variation.
Quote:

dbSNP only merges rs numbers that have an identical set of mappings to a genome and the same allele type (e.g. both must be the same variation type and share one allele in common).  We therefore would not merge a SNP and an indel (insertion/deletion) into a single rs number (different variation classes) since they represent two different types of mutational "events".

In your case, rs6010164 is of the Single Nucleotide Variation (SNV) variation type, while rs267607236 is of the Small deletions or insertions (DELINS) variation type.
Quote:

Despite its name, RefSNP is assigned to all variation types listed below with precise locations for both common and rare variations, including mutations. Most are typically small variations (<= 50bp).

Single nucleotide variation (SNV)
Short multi-nucleotide changes (MNV)
Small deletions or insertions
Small STR repeats
retrotransposable element insertions


text:
igv_plotter error: Xvfb did not start
I am trying to get screenshots from alignments of reads to a genome in IGV (Integrated Genome Viewer) using igv_plotter.
igv_plotter --genome ref.fa readsSorted.bam --locus-file locus_file.txt --max-memory 10G

Traceback (most recent call last):
  File "/Users/user/anaconda/envs/python3env/bin/igv_plotter", line 140, in <module>
    igv_robot.execute()
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/igv_api.py", line 256, in execute
    self._execute_impl(self.__command_queue)
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/igv_api.py", line 387, in _execute_impl
    fake_display.start()
  File "/Users/user/anaconda/envs/python3env/lib/python3.6/site-packages/xvfbwrapper.py", line 79, in start
    raise RuntimeError('Xvfb did not start')
RuntimeError: Xvfb did not start

locus_file.txt only contains this: I:10,000-20,000. From reading online, it seems that this Xvfb is what I need to be able to use IGV without a GUI. I am working on Mac High Sierra and I don't know if I can install Xvfb on Mac OS X or how I can do it. I was quickly out of my depths on the online posts about Xvfb on Mac. I checked that I installed igv_plotter correctly by successfully pulling out the help with igv_plotter -h. I could not find Xvfb using brew.
Alternatively, is there another programmatic way of taking images of aligned reads to a reference, similar to the one below?
Opening IGV and making a screenshot is not a viable option for me, as I would like to export hundreds or thousands or regions.

Hi, What is your question getting igv_plotter to work or find an alternative? If the first, why don't you use the gui to make the screenshots?

@Llopis Maybe because the result is unusable (importing it in Illustrator yields rubbish that’s very time-consuming to unravel).

Did you install XQuartz on your macOS? It should include Xvfb.

@Llopis, I want to automate this for thousands of regions, so opening the GUI is not an option.

@KonradRudolph XQuartz is installed on my macOS. Not sure what it is for but I can open it and it works like a Terminal.


text:
Changing the legend names
ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) + geom_boxplot() +
        scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3")) + theme_classic() +
        labs(title="Comparison",x="Norm", y = "counts") +
        scale_fill_continuous(labels = paste("condition1", "condition2", "condition3"))

I am trying to change the labels of the legend. Currently, it shows the labels that are present in my dataset. However, I want to change it to "condition1", "condition2", "condition3". Above is the code I am trying. It gives me following error.
Scale for 'fill' is already present. Adding another scale for 'fill', which will replace the existing scale. Error: Discrete value supplied to continuous scale In addition: Warning message: Removed 18 rows containing non-finite values (stat_boxplot).

Could someone help me here?

I'm voting to close this question as off-topic because it is not about bioinformatics. Also this problem is answered in Stack Overflow

You don't have any continuous fill scales, since you're filling according to Condition, which is a factor. This is leading to the Error: Discrete value supplied to continuous scale error. You just need to combine your previous scale_fill_manual(values=...) command with the labels= part of your incorrect scale_fill_continuous() command:
ggplot(dataset, aes( x = Norm, fill = Condition, y=log2(measure))) +
    geom_boxplot() +
    scale_fill_manual(values=c("brown1","darkolivegreen4","burlywood3"),
                      labels = paste("condition1", "condition2", "condition3")) +
    theme_classic() +
    labs(title="Comparison",x="Norm", y = "counts")

Something along those lines should work.


text:
Is it possible to use SNP heterozygosity as a proxy for Indel heterozygosity?
I have estimated genome-wide heterozygosity levels using maximum likelihood and classical substitution model implemented in package atlas. These estimates are way more robust than classical SNP calling and simply counting number of heterogeneous locations because they integrate information from a big window into one estimate. I.e. I am quite sure about heterozygosity levels of SNPs of my genomes.
However, now I would like to know also about heterozygosity of indels and call individual snps using my robust estimates as a prior. GATK allows to specify priors on heterozygosity both of SNPs and indels, which might improve variant calls a lot, because default GATK is parametrised for human (Default priors for heterozygous SNPs is 0.001 and for indel 1.25E-4.). But, I work with insects I know that my SNPs heterozygosity is somewhere between 0.01 and 0.03 (1% - 3%).
So here comes the question, if I know the heterozygosity of SNPs, can I simply guess a reasonable prior for the heterozygosity of indels? GATK defaults for priors are 10x different, should I assume the same relation? Is there a reasonable way how to guess a prior for heterozygosity of indels?

SNPs, insertions, and deletions are generated by different biological processes (see for example here). I would expect the per-generation rates of these mutational processes to vary across species. I would also expect the degree to which each mutational process varies across species to vary as well. I think the answer to your question is "probably not without further information".
However, priors are a feature of Bayesian statistics. In Bayesian inference, it is good practice to redo one's analysis using several priors to establish how sensitive one's results are to the chosen prior. I would suggest that you try to call indels using a range of priors and compare the quality of your call sets to determine how important the choice of prior actually is.

1, The processes that generate mutations are indeed different, but the processes maintaining them are the same (drift, selection).
2, Looking at result by different prior is sort of equivalent to using maximum likelihood, no? What would be the point of priors if we had no way how to specify them meaningfully?

I am starting to think, that it's just not a good question because it does not have a clear "yes/no" answer.

Maybe it's worth rephrasing the question to "how do I..." or "What information do I need in order to..."

I thought that posterior probability is likelihood times prior. Then trying multiple priors (or taking a non-informative one) is kind of the same thing as just calculating maximum likelihood estimate without the prior. Or have I missed something?

I see. ML (together with the bootstrap!) and Bayesian inference are both statistical methods that are commonly used to solve similar problems. In that regard, they are "kind of the same thing".  As to the second question, "What would be the point of priors if we had no way how to specify them meaningfully?": Of course, there would be none ;) But the point is that sometimes we can, and the Bayesian approach works regardless of whether you have a strong prior or not.