text:
Lisp in bioinformatics
After many years of programming in Python, C/C++ and other mainstream languages, I am looking for a new language of pleasure, mainly as a part of mental hygiene. Since I really admire the ideas behind Lisp, including the original paper from the 1950's, I would like to learn some Lisp dialect.
Which Lisp dialects are currently suitable for bioinformatics (e.g., for developing simple command-line tools)? I have the following requirements:

Portability across main Unix platforms and reproducibility (this seems to be a huge issue in many Lisp-like languages)
Resulting programs being installable using Bioconda (it should be technically possible at least in future if not now)
Negligible starting time (this disqualifies Clojure, which requires approximately 10 seconds)

I like the idea but I’m not aware of a strong Lisp presence anywhere in bioinformatics. I think a good starting point would be the paper How the strengths of Lisp-family languages facilitate building complex and flexible bioinformatics applications by Khomtchouk &al (and the Hacker News discussion).
Excluding Clojure from consideration will unfortunately close off a great many modern resources (although I fundamentally agree with your technical reason for doing so). Apart from this, Scheme (in particular R7RS) is probably the most-used Lisp dialect, ahead of Common Lisp. Or at least the most actively developed one. As Pierre mentioned, there’s also GNU Guile but, to be honest, despite it being the default for GNU software tools I don’t think it’s frequently used outside GDB.
That said, Clojure’s startup time is not anywhere near 10 seconds on my system. It’s not exactly instantaneous but not much worse than other dynamic languages, despite being backed by the rather sluggish JVM.

Clojure's slow startup time is no longer thanks to Babashka: https://book.babashka.org/

Which LISP dialects are currently suitable for bioinformatics (e.g., for developing simple command-line tools)? I have the following requirements:

I have recently pushed a PR for samtools (mostly for fun) https://github.com/samtools/samtools/pull/752 that uses the GNU Guile library https://www.gnu.org/software/guile/
./samtools view   -g \
   '(define read-filter (lambda (R) (and (> (hts-read-length R ) 10 ) (integer? (string-contains  (hts-read-seq R) "GCCTTGGCTCTTGTTCC")))))' input.bam


text:
Cross-reference with PDB database
I have a list of several thousand proteins and their UNIPROT IDs. I'm looking for an efficient method of cross-referencing it against the PDB tertiary structure database, and get a list of those proteins with a tertiary structure in the PDB database.
I've tried to BLASTP the list of UNIPROT IDs against the PDB database, using the NCBI BLAST portal but encountered too many errors of "Error: Failed to read the Blast query: Sequence ID not found", making the process of manual filtering not convenient and not efficient.
Any ideas?
Thank you!

The UniProt website offers an ID mapping tool that tends to work pretty well:

Then, hit "submit" and you will get the list of corresponding PDB IDs:

Thank you, it worked.
Please aid a PDB-Noob like myself to understand the results: https://ibb.co/jTMGmb
How come each of the UNIPROT ID received several different PDB IDs?
https://ibb.co/b5fRmb

@AdrianSmith that's how PDB works. You'll have one structure for the active site, another for the active site bound by foo, one for the active site bound by bar, another structure for the protein in complex with something else. . . This should be clear enough if you just click on the results you got and see their descriptions.

Thanks for the explanation. I'm looking for a specific protein-target structure (I have a list of ligand-protein interactions and i'm looking for the 3d structure of the binding site if available). Is there any way to search through all of the different "structures"? Does the PDB site offer such a curated database?

@AdrianSmith no idea, I'm afraid. I've almost never used PDB. I suggest you ask a new question explaining what you need in more detail. Ideally, give an example UniProt ID and explain what you would want to see from it.

You can use the SIFTS mapping: http://www.ebi.ac.uk/pdbe/docs/sifts/
From an article:

The Structure Integration with Function, Taxonomy and Sequences
  resource (SIFTS; http://pdbe.org/sifts) is a close collaboration
  between the Protein Data Bank in Europe (PDBe) and UniProt. The two
  teams have developed a semi-automated process for maintaining
  up-to-date cross-reference information to UniProt entries, for all
  protein chains in the PDB entries present in the UniProt database.
  This process is carried out for every weekly PDB release and the
  information is stored in the SIFTS database.

It has no web interface for querying. You just download CSV (or TSV) file with all the mapping data.

Thank you. I'll try it if the ID mapping tool won't provide.

source: Stack Exchange

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

sort: votes


text:
Infer missing UTR features in GFF3 file
I am working on a GFF file that is missing the 5'UTR and 3'UTR information. For example:
  ctg123 . gene      1050  9000  .  +  .  ID=gene00001;Name=EDEN
  ctg123 . mRNA      1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
  ctg123 . exon      1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
  ctg123 . exon      3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
  ctg123 . exon      5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
  ctg123 . exon      7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
  ctg123 . CDS       1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1

Is there a way to add rows in the GFF with the 5'UTR and 3'UTR ranges?

Would we also need to deal with the negative strand? I mean, can you have gff lines referring to the negative strand and, if so, can we assume they will be correct? That the positions will be given with respect to the + strand and the end position will be "before" the start one?

yes, we do have to deal with the negative strand, too. Here is just a simplified example.

Then please [edit] your question and include an example of the negative strand as well. We would also need to know if each gff will have a single entry or if you can have multiple genes in the same gff. The more specific information you give us, the likelier it is that we will be able to give you a useful answer. However, bear in mind that UTRs can be spliced, so this is not a good way of getting the information you need. Getting a full annotation from wherever you downloaded this would be better.

Annoyingly, you can load such a GFF3 file into R (using the GenomicFeatures package) and access the UTRs, but they don't then get saved if you use the export function from rtracklayer (this would have been the simplest solution).
Since I don't know of anything to do this off-hand, I spent a couple minutes to write one. The usage would be addGFF3UTRs.py input.gff output.gff. You can modify a number of the labels to appear as you'd like. This hasn't been extensively tested, but it seems to produce the correct results in a small test I just ran.

The code is hosted on GitHub at the link I gave above and I also include it here:
#!/usr/bin/env python
import argparse

def parseAttributes(kvps):
    d = dict()
    for kvp in kvps:
        k, v = kvp.split("=")
        d[k] = v
    return d

def parseGFF3(fname, topLevel="gene", secondLevel="mRNA", CDS="CDS", exon="exon"):
    f = open(fname)
    genes = dict()
    transcripts = dict()
    g2t = dict()  # A dictionary with gene IDs associated to their transcript IDs

    for line in f:
        if line.startswith("#"):
            continue
        cols = line.strip().split("\t")
        attribs = parseAttributes(cols[8].split(";"))
        ID = attribs["ID"]

        if cols[2] == topLevel:
            genes[ID] = line
            g2t[ID] = []
        elif cols[2] == secondLevel:
            parent = attribs["Parent"]
            g2t[parent].append(ID)

            # Generate a list of [line, [exon entries], [CDS entries]]
            if ID not in transcripts:
                transcripts[ID] = [line, [], []]
            else:
                transcripts[ID][0] = line
        else:
            parent = attribs["Parent"]
            if parent not in transcripts:
                transcripts[parent] = [None, [], []]
            if cols[2] == exon:
                transcripts[parent][1].append((int(cols[3]), int(cols[4]), line))
            elif cols[2] == CDS:
                transcripts[parent][2].append((int(cols[3]), int(cols[4]), line))
    f.close()

    return genes, transcripts, g2t

def findUTRs(transcripts, fivePrime = True):
    for k, t in transcripts.items():
        # If there are no CDS entries then skip
        if len(t[2]) == 0:
            continue
        # Get the strand ("." will be treated as "+")
        strand = t[0].split("\t")[6]

        # For 3' UTR, just swap the strand
        if not fivePrime:
            strand = "+" if strand == "-" else "-"

        exons = [(s, e) for s, e, _ in t[1]]
        CDSs = [(s, e) for s, e, _ in t[2]]
        exons.sort()
        CDSs.sort()
        UTRs = []
        if strand != "-":
            final = CDSs[0][0]
            for s, e in exons:
                if e < final:
                    UTRs.append((s,e))
                elif s < final:
                    UTRs.append((s, final - 1))
                else:
                    break
        else:
            final = CDSs[-1][-1]
            for s, e in exons:
                if e < final:
                    continue
                elif s > final:
                    UTRs.append((s, e))
                else:
                    UTRs.append((final + 1, e))
        t.append(UTRs)

def saveGFF(oname, genes, transcripts, g2t, fivePrime, threePrime):
    o = open(oname, "w")
    o.write("##gff-version 3\n")
    for geneID, geneLine in genes.items():
        o.write(geneLine)
        # Go through each transcript
        for transID in g2t[geneID]:
            t = transcripts[transID]
            if t[0] is None:
                continue
            o.write(t[0])
            cols = t[0].strip().split("\t")

            # Write the exons, then CDS, then 5'UTR, then 3'UTR
            for exon in t[1]:
                o.write(exon[2])
            for CDS in t[2]:
                o.write(CDS[2])
            for idx, UTR in enumerate(t[3]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], fivePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.fivePrimeUTR{};Parent={}\n".format(transID, idx, transID))
            for idx, UTR in enumerate(t[4]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], threePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.threePrimeUTR{};Parent={}\n".format(transID, idx, transID))
    o.close()

parser = argparse.ArgumentParser(description="Parse a GFF3 file lacking UTR entries and add them in. Note that this program assumes a reasonably formatted GFF3 file")
parser.add_argument("--fiveUTRname", default="five_prime_UTR", help="The label for 5' UTR entries (default: %(default)s)")
parser.add_argument("--threeUTRname", default="three_prime_UTR", help="The label for 3' UTR entries (default: %(default)s)")
parser.add_argument("--topLevelID", default="gene", help="The 'type' designating the top-level entry (default: %(default)s)")
parser.add_argument("--secondLevelID", default="mRNA", help="The 'type' designating the second-level entry, typically something like 'mRNA' or 'transcript' (default: %(default)s)")
parser.add_argument("--exonID", default="exon", help="The 'type' designating exons (default: %(default)s)")
parser.add_argument("--CDSID", default="CDS", help="The 'type' designating CDS (default: %(default)s)")
parser.add_argument("input", metavar="input.gff3", help="Input file name")
parser.add_argument("output", metavar="output.gff3", help="Output file name")
args = parser.parse_args()

genes, transcripts, g2t = parseGFF3(args.input, topLevel=args.topLevelID, secondLevel=args.secondLevelID, exon=args.exonID, CDS=args.CDSID)

# Add the UTRs
findUTRs(transcripts)
findUTRs(transcripts, fivePrime=False)

# Write the output
saveGFF(args.output, genes, transcripts, g2t, args.fiveUTRname, args.threeUTRname)

Among other things, the canon-gff3 program in the AEGeAn Toolkit infers and prints missing features like UTRs.
$ canon-gff3 < genes.gff3 > genes-with-utrs.gff3

You can use the add_utrs_to_gff.py script provided by NCBI RefSeq for this purpose. Running the script on your example produces the following:
./add_utrs_to_gff.py example.gff3
ctg123  .  gene             1050  9000  .  +  .  ID=gene00001;Name=EDEN
ctg123  .  mRNA             1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
ctg123  .  exon             1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
ctg123  .  exon             3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
ctg123  .  exon             5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
ctg123  .  exon             7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
ctg123  .  CDS              1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  five_prime_UTR   1050  1200  .  +  .  ID=utr00002;Parent=mRNA00001
ctg123  .  three_prime_UTR  5301  5500  .  +  .  ID=utr00004;Parent=mRNA00001
ctg123  .  three_prime_UTR  7000  9000  .  +  .  ID=utr00005;Parent=mRNA00001

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

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


text:
"No records to view" error in the "Drug and Drug Target Mapping" PDB's search engine
I'm using the "Drug and Drug Target Mapping" search engine (http://www.rcsb.org/pdb/ligand/drugMapping.do).
I've noticed that there are some drugs (generic names) that if they appear in the query, the search results end in "No records to view".
I've failed to find a pattern in identifying those drugs.

The site is no longer active.
One replacement is the 'anatomical therapeutic chemical' (ATC) classification. Information is found on RCSB, PDB,
https://www.rcsb.org/docs/search-and-browse/browse-options/atc
...and WHO, https://www.whocc.no/atc/structure_and_principles/
The question remains pertinent.

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


text:
Can Mauve export indels or is there an whole-genome aligner that can export indels?
I have aligned two whole genomes of the same species in Mauve using progressiveMauve. One of the genomes is 2 megabases longer than the other so I know it must contain some insertions. Looking at the alignment in the Mauve GUI, it seems that insertions are small and scattered throughout the genome. I have exported the SNPs from the GUI. Is there a similar option for indels?
I have also inspected the Mauve file written after the alignment. It contains the two original sequences in FASTA format but I noticed in contains gaps denoted by a dash (-). Would it be possible to create a VCF from that file?


text:
how to count the number of each allele of CNVs?
I have a bunch of files (vcf format) each containing site and sample annotations of a copy number variant (following example is part of one of my files, I am showing 11 individuals but they are nearly 900 individuals)
16  32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3

I want to count each allele of CNV in each individual(sample) if there are 0, 1 or 2 copies!
this is my desired output (rows are individuals):
CN0 CN2 CN3 CN4 CN5
0   0   2   0   0
1   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   2   0   0
0   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   0   0   0
0   0   2   0   0
1   0   1   0   0

Any idea how to do that?

Cold you explain your file format a bit? How do we get from the VCF line you show to the table of your desired output?

just edited my post to clarify that it is a vcf form file. I used vcftools --poition option to extract the position that I needed! and how to create the table is exactly the question that I have .....!

Here's a pure python solution to process your example. It should be straightforward to loop over all the rows of your vcf files
my_line='16 32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3 '
my_line=my_line.strip() # remove trailing whitespace
my_fields=my_line.split()
names=my_fields[4].split(',') # extract CNV names
names=[x[1:-1] for x in names] # remove '<>' symbols
print('\t'.join(names))
for sample in my_fields[9:]:
    genotypes=[0]*(len(names)+1) # initialize array of zero counts, including reference
    for allele in sample.split('|'):
        genotypes[int(allele)]+=1
    print('\t'.join([str(x) for x in genotypes[1:]])) # print all columns except reference

Thanks @heathobrien, this may work but in the example above I just showed genotypes for a few individuals, I have nearly 900 individuals and It is not convenient to copy all genotypes in your code in "my_line" part! Do you have any alternative to be able to fix this?

You'll probably want to include the code I wrote as part of a larger script that reads your genotypes file and writes the output to another file rather than printing it to the screen. This isn't a script writing service, so take a look at an python I/O tutorial or the documentation for the csv module (or pysam)


text:
Where can I find a list of SNP rs IDs that belong to the X chromosome?
Where can I find a list of SNP rs IDs that belong to the X chromosome?
I'm trying to exclude the SNPs on X-chromosome in a data set, but I don't know where these SNPs come from.
I tried to use dbSNP to search, it returns 26933249 SNPs, but the webpage seems to refuse me to download the information.

Have you tried the send to option in the top right corner, just below the search bar?

In what species? Your search didn't specify a species, are you really interested in rsIDs from any species?

You just searched for X[Chromosome], you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for X[Chromosome] human[organism]. Then, click on the send to link at the top left and download the data:

This is veeeeeery slow, however. I suggest that you instead download the data for chrX from dbSNP's FTP site (this assumes your species of interest is human): ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606/chr_rpts/chr_X.txt.gz
Once you have this file, you can use it to filter your dataset:

Collect all rsIDs from the file:
zcat chr_X.txt.gz | awk '/^[0-9]/{print "rs"$1}' > chrX.rsIDs

Filter your dataset.
You haven't shown the format of your dataset, so if it is just a list of rsIDs like this:
$ cat dataset
rs77669712
rs199706910
rs2799940
rs12205639
rs10958205
rs200778574

You can remove any that are on chrX with:
$ grep -vxFf chrX.rsIDs dataset
rs77669712
rs199706910
rs12205639
rs10958205
rs200778574

That will print all rsIDs from the text file dataset that are not also in the text file chrX.rsIDs. Note the -x option for grep which tells it to only consider matches that cover the entire line. This ensures that rs1234 is not removed when you're looking for rs123.


text:
bed file to gtf conversion using awk
I have a bed file like this:
chr1    12227   12612   LOC102725121    0   +
chr1    12721   13220   LOC102725121    0   +
chr1    12227   12612   DDX11L1 0   +
chr1    12721   13220   DDX11L1 0   +
chr1    14829   14969   WASH7P  0   -
chr1    15038   15795   WASH7P  0   -
chr1    15947   16606   WASH7P  0   -
chr1    16765   16857   WASH7P  0   -
chr1    17055   17232   WASH7P  0   -
chr1    17368   17605   WASH7P  0   -

I am trying to convert into gtf file for downstream analysis. Gtf file should look like this:
chr1   myIntrons   feature 3073253   3074322   .   +   .     gene_id "intron1"; transcript_id "intron1"
chr1   myIntrons   feature 3075858   3079636   .   +   .     gene_id "intron2"; transcript_id "intron2"

For the gene_id part in the 9th column in the double quotes I want my 4th column and in transcript_id I want combination of column1,column2,and column3 so that my output looks like this:
Showing for the first line
chr1    myIntrons   intron  12228   12612   .   +   .   gene_id "LOC102725121"; transcript_id "chr1:12227-12612"

I am using this code but under double quotes it is not expanding $1 $2 etc
awk -v OFS="\t" '{print $1,"myIntrons","intron", $2+1, $3, ".", $6, ".", "gene_id $4; transcript_id $1:$2-$3"}' file.bed| head

Instead of the "gene_id ..." bit, you want sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"", $4, $1, $2, $3).


text:
Is there a tool that can take a protein's amino acid sequence and would display it's locus on the genome?
I have the UNIPROT IDs, PDB IDs and FASTA files of several known proteins. I am looking for a tool that can take as input the protein's amino acid sequence and display the coding nucleotides of those amino acids. Moreover, I look for a tool that would be able to display those nucleotides' location on the human genome.
The final outcome I desire to achieve is to find SNPs in the exons of the protein.
A tool that can handle several amino acid sequences at once is preferred for I have plenty of proteins. I've tried to use UCSC's BLAT human genome search but it was too cumbersome. I've also tried BLAST, but I'd like a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.
Thank you for your assistance.

For something "straightforward" you'd need the correct reading frame, strand orientation, and splicing to get back to a nucleotide sequence to align to the genome. The tools mentioned should help take a lot of work out of this.

I think the mistake here is thinking that you need to use the sequence in order to find the genome co-ordinates. If these proteins come from a well annotated speices, you can easily do this with BioMart.

Go the the ENSEMBL web site.
From the top menu, select BioMart
Under "-CHOSE DATABASE-" select Ensembl Genes, and under "-CHOOSE DATASET-" select your species (eg. Human Genes).
Select "Filters" from the lefthand menu and then expand the "GENE" section.
Select the "Input external references ID list [Max 500 advised]" and from the drop-down select "PDB IDs" or the correct "Uniprot KB" ID.
Paste your IDs into the box.
Select "Attributes" on the left-hand menu. Making sure you've the "Fetaures" radio box select, you can select Chromosome and either Gene Start/Gene End or Transcript start/Transcript end from the GENE section. These coordinates will also contain the untranslated regions unfortunately.  You may also want to play about with the various things under "PROTEIN DOMAINS AND FAMILIES", particurlarly under "Protein Features", although I'm not sure what they all mean.
When you are happy with the selection, select the big "Results" at the top to get the coords.

You can also get out Sequences for both the full transcript and the CDS and known germline and somatic variants that fall within the boundaries of the genes.

Thank you, very useful tool. Can it display the external results references IDs (ie the PDB IDs)  at the results page?

Have you tried BLAST, or specifically tBLASTn?  This will take a protein and compare it any set of reference genomes.  You can also BLAST multiple sequences at once.

I've tried it, but I look for a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.

I think that you will have to define "straightforward" if you want help.  BLAST and BLAT are the standard tools.  You need to describe what issues you have with those tools and what you would like instead.


text:
How to concatenate "by chromosome"-VCFs?
I have a several VCFs which are VCF which only contain information by chromosome. That is, there's a chromosome 1 VCF (with only chr1), a chromosome 2 VCF (with only chr2), etc.
I checked to make sure that these VCFs were valid via VCFtools, i.e.
$ vcf-validator chr1.vcf

which works---these are valid VCFs I was given.
Now, I would like to combine these VCFs into one VCF.
I naïvely tried the following cat operation:
$ cat chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf > total_chroms.vcf

This doesn't work properly though
$ vcf-validator total_chroms.vcf
The header tag 'contig' not present for CHROM=chr1. (Not required but highly recommended.)
Could not parse the line, wrong number of columns: [##fileformat=VCFv4.2\n]
 at /path/vcftools-0.1.14/perl/Vcf.pm line 172, <__ANONIO__> line 191016.
     Vcf::throw('Vcf4_2=HASH(0x1ae7208)', 'Could not parse the line, wrong number of columns: [##filefor...') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 335
     Vcf::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 3457
     Vcf4_1::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 2574
     VcfReader::run_validation('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 60
     main::do_validation('HASH(0x16ada68)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 14
$

What tools are available to merge these VCF's together into a total VCF?

I would recommend bcftools concat.  You can't just cat them together because each file has a header section.  The bcftools command will handle all that for you.  Each vcf file must be sorted prior to calling concat
bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf

It's a bit odd. vcf-validator chr1.vcf shows no problems. I believe the output of bcftools concat -o chroms12.vcf chr1.vcf chr2.vcf shows an error: FIXME: sequence name chr1 in chr1.vcf

@ShanZhengYang the errors seem to indicate that you have chr1 in the CHROM field for some rows, but no corresponding line in the header. Can you check if that is true in the chr1.vcf file? Maybe the validator somehow is missing this error?

@juod The VCF header looks correct to me. #CHROM POS ID REF ALT QUAL FILTER INFO

@juod "but no corresponding line in the header" Maybe I don't understand this. Besides the code pasted above, there is no other "chrom" information in the header

According to the VCF spec https://samtools.github.io/hts-specs/VCFv4.1.pdf there should be header lines for each chromosome this is from the first example in the linked documents - ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>

@Bioathlete I see. I think I understand now---this issue has been discussed before: https://github.com/samtools/bcftools/issues/326
I know that these VCFs are from hg38. Could I simply add the line ##contig=<ID=chr1,length=248956422>, etc. to the VCFs in question? I would just take the chromosome lengths from hg38

@ShanZhengYang try that - it should work. I don't even think it's necessary to provide the length, but you could experiment and let us know.

@juod The above worked! Thanks for the help!

The best tool for the job is probably bcftools as suggested by Bioathlete, but you can also do this manually. You just need to collect all header lines from all vcf files, remove any duplicates, and then print all the headers + the actual data into the new one:
grep '^##' chr*vcf | sort | uniq > all.vcf
grep -m1 '^#CHR' chr1.vcf >> all.vcf  ## Get the chr header line
grep -v '^#' chr*vcf >> all.vcf

I appear to get a malformed VCF header with the above commands

@ShanZhengYang malformed how? Could you edit your question and show us an example of your files?

So, I see the above grep commands give the same behavior with all of my standard VCFs. I'll get something like this in all.vcf:
chr10.vcf:##contig=<ID=chr10,length=133797422> chr10.vcf:##fileDate=20121011 chr10.vcf:##fileformat=VCFv4.2 chr10.vcf:##reference=path/reference.fa chr10.vcf:##source=foobar chr11.vcf:##contig=<ID=chr10,length=133797422> chr11.vcf:##fileDate=20121011 chr11.vcf:##fileformat=VCFv4.2 chr11.vcf:##reference=path/reference.fa chr11.vcf:##source=foobar ....

use picard GatherVcf: http://broadinstitute.github.io/picard/command-line-overview.html#GatherVcfs

Gathers multiple VCF files from a scatter operation into a single VCF file. Input files must be supplied in genomic order and must not have events at overlapping positions.

The error I get here is Exception in thread "main" picard.PicardException: In order to index the resulting VCF input VCFs must contain ##contig lines.. It is true, my "by chromosome"-VCFs do not contain a ##contig line in the header...but I'm not sure that would matter

Edit above: Apparently this issue has been discussed before: github.com/samtools/bcftools/issues/326
Question: how could I add these contig lines if I don't know this information?I know that these VCFs are from hg38. Could I simply add the line ##contig=<ID=chr1,length=248956422>, etc. to the VCFs in question? I would just take the chromosome lengths from hg38

"but I'm not sure that would matter" it checks the order of the contigs using the dict. use UpdateVcfSequenceDictionary https://broadinstitute.github.io/picard/command-line-overview.html

"it checks the order of the contigs using the dict. use UpdateVcfSequenceDictionary". Ah, I understand now. Thanks for the help!


text:
How can I count the number of reads that support a variant in a bam file?
I am calling variants from a human sample using bwa mem to align the reads and gatk to call the variants. I'm trying to understand why a specific variant was not called in my sample. I have checked the bam alignments in a GUI viewer and I can see that there are reads supporting the missing variant. It looks like the issue is a low allelic balance, with far more reads supporting the reference than the alternate allele but I want to get the actual numbers.
So, given a specific variant like this:
chr22   425236  C T

How can I count the number of reads in my sample.bam file that support that variant and the number that don't on Linux?

If you don't mind a bit of manual counting, then samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam will produce output where you can count the bases for that position. You could, of course, use the command line to do most of that automatically:
samtools mpileup -f reference.fa -r chr22:425236-425236 alignments.bam |
    cut -f 5 | tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c

That'll give you a count of how many of each base were seen.

ASCIIGenome (I'm the author) has a command, filterVariantReads, designed to inspect reads having a variant at a position or range. It would go along these lines:
ASCIIGenome -fa genome.fa aln.bam

Then go to the region of interest and use:
goto chr9:4917981-4918161
filterVariantReads -r 4918011

From this:

You get:

You can also print to screen or save to file the underlying sam records with:
print > reads.sam

For several regions, the whole process can be scripted for automation.
Hope this helps!

Hello, I am having a trouble to download this tool neither with conda nor from the source code. The .jar file is missing...

@user3224522 What have you tried? The zip file ASCIIGenome-1.16.0.zip should contain what you need. conda install asciigenome should also work. Please submit an issue to GitHub or post a question with more detail.

The variant record in the VCF produced by GATK should include information about depth -- both total read depth (DP) and allelic depth (AD), ie depth per allele.
Keep in mind however that the caller makes its decision based on more than just allelic depth. There are some additional statistics involved, including the quality of the individual bases. Also, if any filtering was applied, your variant may have been excluded due to other statistics, like for example if there is evidence of strand bias -- eg if all the reads supporting your variant are in the same orientation, that is a classic red flag suggesting it's an artifact.


text:
bed entries collapse based on score column
I have a bed file where some entries have the exact same chr:start-end but change in the name and score column. I would like to apply bedtools or a similar tool to obtain a filtered file where for each set of chr:start-end entries, only the one with the highest score (column 5) is used, the others filtered out. I would like to do this with a bed format toolkit (bedtools, etc) rather than scripting it.

I searched for bedtools merge or bedtools collapse, but couldn't find a way.

Use BEDOPS bedmap --max-element or bedmap --max-element-rand to grab the maximal-scoring element over mapped intervals, e.g.:
$ sort-bed elements.unsorted.bed5 > elements.bed5
$ bedops --merge elements.bed5 > merged.bed3
$ bedmap --max-element merged.bed3 elements.bed5 > answer.bed

The file answer.bed will contain the maximal-scored element from elements.bed5 that overlaps the merged intervals in merged.bed3.
The --max-element and --max-element-rand options differ on handling ties. The first option grabs the first match from ties by sort-bed (lexiographical) sort order; the second option grabs a random match on ties.
The option you pick would depend on what you're doing with the answer. Picking a random hit might be better if you're doing downstream statistical tests, for instance.
You can add --sci and --prec <int> options to format the score output to your desired level of precision. I think the default assumption is six digits of precision; adjust as needed.
There are also --min-element and --min-element-rand options for, conversely, reporting minimal-scoring elements.
For more information: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#score-operations

groupby from bedtools might help you:
bedtools groupby -i variantsToRepeats.bed -grp 1-3 -c 5 -o max

-grp to select the columns to group by
-c to select the columns to summarize
-o to select how to summarize the column selected.

Could you please expand your answer? An explanation might help future readers. Thanks! I edited your question from the help page you linked. I hope it is now better, and I understood it correctly

How lazy can people be? I even linked the help page and to anybody with an IQ over 80 the options should be immediately obvious.

Without the explanation your answer becomes pretty unhelpful if the location of the help files changes


text:
Make ipyrad use cuda-enabled NVidia card on Ubuntu
I want to use ipyrad on a new Ubuntu machine that has an NVidia Quadro K2000 card with 384 cores. One can configure ipyrad to run on a linux cluster.   Do I have any options to get ipyrad to access these cores?
I'm assuming that I will be able to install cuda drivers on the machine. If ipyrad can't be made to do this, is there another way to do de novo assembly of GBS data that will take advantage of this resource?

GPU core ≠ CPU core. The two aren’t interchangeable.

You would have to get CUDA Python set up on your system and likely modify ipyrad to actually use it (for example, by incorporating numba). So while this is possible, it will require significant effort on your part. In general, there are very few bioinformatics tools that support GPUs.


text:
What is the --rev_comp_mapping_barcode parameter in the QIIME1 script extract_barcodes.py?
I'm using the QIIME1 scripts extract_barcodes.py to extract barcodes from a dual-indexed MiSeq library, and then demultiplex my reads using split_libraries_fastq.py.  I have been having trouble correctly using a couple of parameters relating to reverse complementing barcodes.
When extracting the barcodes with extract_barcodes.py, I use --rev_comp_bc1 (as directed by the individual in charge of sequencing). Then when I demultiplex the reads with split_libraries_fastq.py, sometimes I need to include the parameter --rev_comp_mapping_barcodes to get the script to work, and sometimes I don't.  I'm not sure why, unless the barcodes I'm receiving in my mapping files are sometimes written in different orientations (I don't write them).
I think the parameter --rev_comp_bc1 is self-explanatory, but could someone explain what --rev_comp_mapping_barcodes does and when to use it?
To make things simple we can assume that the two barcodes are written from 5' to 3' in the mapping file.

I dowvoted this question despite being well write because I don't see any reference to what have you done to understand the differences between the parameters. Did you look up to the documentation? Have you asked the responsible of the sequencing why sometimes is like this? When you update the question I will remove the downvote.

Of course I have. There's no point in being critical. You obviously don't know the answer so why bother.

Sorry I was so critical. I did some effort to answer your question and then realized you didn't disclose if you where aware of the documentation and how it hadn't helped you. I was a bit disappointed to have spend some time for nothing, that's why I decided to leave the comment to prevent other to invest more time to answer your question in pages and papers you are already aware. I hope you find the answer somehow. Best wishes

As the help page says:

--rev_comp_mapping_barcodes
      Reverse complement barcode in mapping before lookup (useful if barcodes in mapping file are reverse complements of golay codes) [default: False]

As the help page of extract_barcode says:

--rev_comp_bc1
  Reverse complement barcode 1 before writing [default: False]

You seem to be doing two complete different things. extract_barcode.py is designed to

format fastq sequence and barcode data so they are compatible with split_libraries_fastq.py

While the first one is for:

demultiplexing of Fastq sequence data where barcodes and sequences are contained in two separate fastq files

So it seems that depending on the machine used or the previous processing you are given with the barcode in another file or together. I would advice to ask the person who gave you this data.


text:
Finding the nearest gene to a specified gene region
I've got some data on genomic regions of interest in a bedgraph-like format from the Mus musculus [mm10] genome, and would like to find out the nearest gene to these regions.
Some of the regions will reside within a gene (in which case I'd like to know that they're within a gene), some regions may reside within multiple genes (in which case I'd like to know all gene names), and some regions are not within a gene (in which case I'd like to know the closest gene: the shortest sequence gap between either end of the gene and either end of the tagged genomic region. In this case I want to do an unstranded search, because it's possible that I might have regions describing suppressing transcripts. Here's some example data:
Ref         Start             End      fwd     rev
chrX     73716152        73716152        0      -2
chrX    167207094       167209162        0      -3
chrY     30770844        30772724        3       0
chr1     24613189        24613641        0       3
chr1     24613971        24614812        0       4
chr1     35788729        35789387       -5       0
chr1     93151586        93160225        0      -2
chr1    107511454       107525597       -2       0
chr1    130882067       130887413        0      -3
chr1    134182466       134189942        2       0
chr1    135818954       135832956       -2       0
chr1    149520498       149521157        0      -3
chr1    166166500       166166500       -2       0
chr2      5379059         5379374        0      -2
chr2     32379110        32381837        0      -2
chr2     32381854        32381858        0      -2
chr2     32384633        32384633        0      -2
chr2     32384636        32387773        0      -3
chr2    120539390       120563799        0      -2
chr2    125946785       125975527        2       0
chr2    127633311       127633684        0      -2
chr2    127634841       127656210        0      -3
chr2    164562717       164562717       -2       0
chr2    164562719       164562720       -2       0
chr2    164568498       164568499       -2       0
chr2    181188007       181188388       -2       0

I'd like to add another 'Gene' column to this table that gives the official names of the nearest gene(s) for each region defined by <Ref>:<Start>-<End>.

What would be the desired output format? Would you like the name of the gene in that new column? Are you interested in the nearest gene in the same strand or it doesn't matter? What happen if multiple genes are in those regions?

You're looking for bedtools closest, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.
  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

You'll need to strip the header and ensure the file is sorted. I've piped the output to awk to reformat things. The output then is something like:
chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.

With the prompting from @Devon_Ryan's answer, I came up with a solution that worked well for me. I downloaded the mm10 RefSeq gene locations from UCSC [Fields: chrom, txStart, txEnd, name2], and sorted both input files using bedtools sort. From there the bedtools closest command got me almost to where I wanted to go:
(echo -e "Chr\tStart\tEnd\tfwd\trev\tm_Chr\tm_Start\tm_End\tGene\tDistance";
  bedtools closest -d -a regions.bed -b NCBI_RefSeq_mm10.bed) > annotated.tsv

Output example:
Chr     Start   End     fwd     rev     m_Chr   m_Start m_End   Gene    Distance
chrX    73716152        73716152        0       -2      chrX    73686177        73716204        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73717815        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716401        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716502        Bcap31  0
chrX    167207094       167209162       0       -3      chrX    167207092       167213479       Tmsb4x  0
chrX    167207094       167209162       0       -3      chrX    167207093       167209218       Tmsb4x  0
chrY    30770844        30772724        3       0       chrY    30725674        30738546        Gm29582 32299
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    35788729        35789387        -5      0       chr1    35869594        35881164        1110002O04Rik   80208
chr1    93151586        93160225        0       -2      chr1    93151353        93160935        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151353        93160957        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160948        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160870        2310007B03Rik   0
chr1    107511454       107525597       -2      0       chr1    107500891       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511422       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511509       107525600       Serpinb2        0
chr1    130882067       130887413       0       -3      chr1    130882073       130887411       Il24    0
chr1    134182466       134189942       2       0       chr1    134182169       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182403       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182408       134190031       Chil1   0
chr1    135818954       135832956       -2      0       chr1    135818597       135833341       Lad1    0
chr1    149520498       149521157       0       -3      chr1    149410728       149449931       Gm29398 70568
chr1    166166500       166166500       -2      0       chr1    166130459       166166510       Gpa33   0
chr2    5379059         5379374         0       -2      chr2    5293456         5676046         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5383533         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5569754         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5601302         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5715350         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5298992         5715357         Camk1d  0
chr2    32379110        32381837        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32381854        32381858        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32384633        32384633        0       -2      chr2    32384636        32387739        Lcn2    3
chr2    32384636        32387773        0       -3      chr2    32384636        32387739        Lcn2    0
chr2    120539390       120563799       0       -2      chr2    120506829       120563831       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563851       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563853       Zfp106  0
chr2    125946785       125975527       2       0       chr2    125859108       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125859228       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866111       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866218       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866219       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866222       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125920565       125983587       Galk2   0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    164562717       164562717       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562717       164562717       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    181188007       181188388       -2      0       chr2    181187221       181188506       Ppdpf   0
chr2    181188007       181188388       -2      0       chr2    181187342       181188504       Ppdpf   0

After that it was just a matter of aggregating the gene names for each unique region with a quick R script:
library(dplyr);
data.tbl <-
    read.delim("annotated.tsv") %>%
    group_by(Chr, Start, End, fwd, rev) %>%
    summarise(Genes = paste(unique(Gene), collapse=","), Distance=min(Distance));

write.table(data.tbl, "collapsed_annotated.tsv",
            sep="\t", row.names=FALSE, quote=FALSE);

Collapsed output example:
Chr     Start   End     fwd     rev     Genes   Distance
chr1    24613189        24613641        0       3       Col19a1 25655
chr1    24613971        24614812        0       4       Col19a1 26437
chr1    35788729        35789387        -5      0       1110002O04Rik   80208
chr1    93151586        93160225        0       -2      2310007B03Rik   0
chr1    107511454       107525597       -2      0       Serpinb2        0
chr1    130882067       130887413       0       -3      Il24    0
chr1    134182466       134189942       2       0       Chil1   0
chr1    135818954       135832956       -2      0       Lad1    0
chr1    149520498       149521157       0       -3      Gm29398 70568
chr1    166166500       166166500       -2      0       Gpa33   0
chr2    5379059         5379374         0       -2      Camk1d  0
chr2    32379110        32381837        0       -2      1110008P14Rik   0
chr2    32381854        32381858        0       -2      1110008P14Rik   0
chr2    32384633        32384633        0       -2      Lcn2    3
chr2    32384636        32387773        0       -3      Lcn2    0
chr2    120539390       120563799       0       -2      Zfp106  0
chr2    125946785       125975527       2       0       Galk2   0
chr2    127633311       127633684       0       -2      Mal     0
chr2    127634841       127656210       0       -3      Mal     0
chr2    164562717       164562717       -2      0       Wfdc2   0
chr2    164562719       164562720       -2      0       Wfdc2   0
chr2    164568498       164568499       -2      0       Wfdc2   0
chr2    181188007       181188388       -2      0       Ppdpf   0
chrX    73716152        73716152        0       -2      Bcap31  0
chrX    167207094       167209162       0       -3      Tmsb4x  0
chrY    30770844        30772724        3       0       Gm29582 32299

Another option is BEDOPS closest-features, which efficiently locates nearest annotations, overlapping or otherwise.
$ closest-features --closest --delim '\t' intervals.bed genes.bed > answer.bed

It's easy to reformat output by piping to awk, say.
If files are very large, you can parallelize work into "map-reduce" form by specifying the chromosome name as an option:
$ closest-features --chrom chr1 --closest --delim '\t' intervals.bed genes.bed > answer.chr1.bed
$ closest-features --chrom chr2 --closest --delim '\t' intervals.bed genes.bed > answer.chr2.bed
...

To reduce, just bedops -u all the answer.*.bed files.
Running these per-chromosome jobs in parallel makes the work take about as long as searching chr1, which tends to be the largest chromosome.


text:
predict the foldability of single-stranded DNA molecules
I have a list of regions of the human genome and I want to predict if single-stranded molecules in a buffer would tend to fold and create pin structures by sequence self-complementarity. What's the most precise software/parameter set I can use to such effect?
TRIED SO FAR:
Vienna RNALFold (unsure if applying an RNA software to DNA is adequate):
Compute locally stable RNA secondary structure with a maximal base pair span.
For a sequence of length n and a base pair span of L the algorithm uses only
O(n+L*L) memory and O(n*L*L) CPU time. Thus it is practical to "scan" very
large genomes for short RNA structures.
Output consists of a list of secondary structure components of size <= L, one
entry per line. Each output line contains the predicted local structure its
energy in kcal/mol and the starting position of the local structure.

Example:
samtools faidx $ref $region | seqtk seq -A - | $HOME/viennarna/bin/RNALfold --noconv
>chr1:100000-100150
.((((.(......).)))). ( -2.70)  127
.(((.((........)).))). ( -5.70)  123
.(((((((.((........)).)))..)))). ( -6.30)  119
.(((.(.(((.((.....))))).)..))). ( -7.00)  117
.((((.((........)).)))). ( -8.30)  109
.((((..((.((..(..(((.....)))..)..)))))))). (-10.10)  103
.((((.((..((..(((...((((.((........)).)))).....)))))..)))))). (-10.50)   90
.((.((.((((....................)))).)).)). ( -4.10)   77
.((((.((.((((....................)))).)).)))). ( -5.70)   75
.((((..((((.((....)).)))).)))). ( -2.40)   72
.((.((..((((..((((.((....)).)))).))))..)))). ( -5.50)   65
.(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).). (-16.20)   40
.((((......(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).).....)))). (-18.10)   30
.(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......))). (-18.50)   23
.((....(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......)))....)). (-18.70)   17
.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))). ( -9.60)    7
.((.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))).)). (-10.90)    4
.(((..(((((........)))))..))). ( -3.60)    1
ACTTAAGTTGTAGAAGGGAATAACGCAAGAGTGAATTTAGGGCGGGGCAAAAGGATAAATTTTACGGTACAAAGTTTCTACGGGTTTATATATGTATAACTAAGTCCAAGCGCGGGGGATATGGCCAGTGCACAACGGCGGGCATCATAAT
 (-21.70)

I think this is a nice question but a bit open-ended. Are you searching for other tools because this didn't convince you ? Or what requirements are you looking in a other tools? How do you define "best"? Less FDR or more precision ?

I edited my question to "most precise" and added that the reason I am unsure about RNALFold is because it's an RNA software applied here to DNA molecules.

Great edit, I didn't realize you were using DNA as input, I thought you were using RNA from those regions

You could also take a look at Infernal.
I think your best bet is to use RNA-focused programs; I'm not aware of a DNA-specific one, which makes some sense because DNA is typically double-stranded, so the secondary structure community is primarily focused on RNA. I would predict that secondary structure potential of single-stranded DNA would be similar to that of RNA.

Thanks for recommending. I learned about Infernal about 5-8 years ago, and after reading the documentation as it is now, it seems to me that there isn't an easy way of asking for the folding of DNA sequences. The commands seem to be based on building alignments. Am I wrong?


text:
Help with generating network from gene expression data
I'm struggling to figure out how to go from gene expression data to  be able to create a network out of the most expressed genes in a cell line. Ok, I have developed an application for Cytoscape where I want to compare my analysis with research on the NCI-60 cell lines. My goal is to do these steps:
1: Use one of the cell lines in the NCI-60 cell lines, in my case I will use the NCI_H23 where there is 3 samples in this dataset.
2: Use the most expressed genes in the NCI_H23 cell line to create a subnetwork using a protein interaction network database like iRefIndex, etc.
After these two steps I should finally be able to upload the network to Cytoscape and perform my analysis, but being a programmer  I'm kind of stuck. I have tried using one of the 3 samples where I mapped the Affymatrix IDs to gene symbol and then selecting the top 20% genes with the highest values using a cutoff value, but I'm not sure if this is the correct way to do it and I'm not sure how to do step nr.2 above if my list of genes is correct.
My application investigates drug combinations based on graph theory so I want to create a network that represents a cell-line so I can compare my scoring of combinations with the ones from the NCI Almanac study. I don't want to compare it with other cell-lines, but rather find the most expressed genes in that cell-line to then retrieve interaction data from iRefIndex to get edges for the network. I guess I should average the 3 samples first before I repeat what I listed above.
I would really appreciate help on this as I've been stuck on this for over a week now.
What I'm stuck at now is trying to create a network from iRefIndex (using iRefR) using the set of gene symbols that I have extracted from the 3 samples. To be honest, I'm finding it hard to find the right tools and assessing what kind of interaction data i should use from iRefIndex etc, so I guess I am looking for confirmation that what I'm trying makes sense.

You seem to make quite a lot of guesses, maybe you should seek guidance with an advisor (offline). Usually one selects the genes that have a change in the expression when comparing with and without the drug. Have you read how other people investigate drug combinations? I don't know why have you selected the iRefIndex database, but are you looking for a interaction database or a pathway database?


text:
Mapping Information from SAM/BAM file
I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM)

How many PacBio reads are mapped to at least one illumina reads
A Table showing each pacbio read, the number of illumina reads that
mapped, and total coverage (i.e. percentage of mapped regions)

By "percentage of mapped regions", do you mean the percentage of the PacBio read that was covered/mapped by illumina reads?

If you sort and index your BAM file:

samtools idxstats. The third column is the number of Illumina reads mapped to each PacBio read. You can then pipe that to awk: samtools idxstats foo.bam | awk '{if($1 != "*" && $3 > 0) tot += 1}END{print tot}. That will print the number of PacBio reads with at least 1 read mapped to them.
You'll want samtools depth for that

samtools depth -aa foo.bam | awk '{if($1 != lread) {if(lread != "") {print lread, 100*lcov/llen}; lread = $1; lcov=0; llen=0}; llen += 1; if($3>0) {lcov +=1}}END{print lread, 100*lcov/llen}'

That won't give you the number of illumina reads per pacbio read, but you can get that from samtools idxstats (column 3) and join the results.
BTW, your usage of "coverage" is non-standard. You want "percentage coverage" rather than "coverage", which actually means "fold-coverage" (i.e., number of reads times their average length divided by the total genome/chromosome/etc. length).
If you instead want "fold-coverage" then the easiest way is to first know the length of your reads. Presuming you have 100 base reads (for paired-end reads, this would be the same, presuming the ends don't overlap):
samtools idxstats foo.bam | awk '{print $1, $3, 100*$3/$2}'

That will give you a table of PacBio read name, number of aligned reads, fold-coverage.

This infomation can be obtained with
samtools index mybam.bam
samtools idxstats mybam.bam

You'll get a table with one row per PacBio read, the length, the number of mapped reads aligned to it and the number of unmapped reads aligned to it.
The number of PacBio reads mapped to at least one illumina read is the number of rows where the 3rd column is greater than 0
samtools idxstats mybam.bam | awk '$3>0' | wc -l


text:
Removing repeated reads from nanopore 1D² reads
This is summary of discussions and a question that was posted by Devon O'Rourke on Twitter
Following Albacore basecalls on a 1D² library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D² script. I'd like to use the reads generated from these scripts in a genome assembly.
After a little bit of grep searching, it seems like there are overlapping reads among the 1D and 1D² directories - namely, those 1D reads that generated 1D² reads!
How do people approach the choice of selecting which directory to pull the unique read from if you're planning on assembling a genome? I'd think that the 1D² reads are of higher quality and be preferred, but perhaps there are 1D² reads which are lower in quality than one of the two 1D pairs?
Just looking for advice on how folks approach 1D² libraries, fastq outputs, and subsequent mapping. For what it's worth I'm hoping to use Minimap2 to proceed with alignments but am happy to hear others thoughts on alternatives.
Update
I tried the first solution proposed by @gringer, but my filtered_1d_reads.fq file contains same number of reads as raw_1d_reads.fq file.
I'm pretty sure the problem is the naming convention difference between the linear and paired reads. linear reads always follow a read_id of alphanumeric (8)-(4)-(4)-(4)-(12), while paired reads are (8)-(4)-(4)-(4)-(20)-(4)-(4)-(4)-(12). When I looked at it in more detail, the paired read "read_id" naming is just a concatenation of the two linear "read_id"'s that are partnered. So that (20) value is really the last (12) of the first linear read, and the first (8) of the second linear pair.

Bit silly, but here it goes. > I'd think that the 1D² reads are of higher quality and be preferred,

but perhaps there are 1D² reads which are lower in quality than one of
the two 1D pairs? Based on this it would make the most sense to have a script preferentially select 1D2 reads unless they have a lower quality than the 1D 'parents', and keep all 1D reads that did not contribute to a 1D2 read. I don't have a dataset to play with, but this shouldn't be too hard to write I believe.

An additional QC step (as suggested by @Wouter_De_Coster) would be a good idea. It could probably be done by using the run statistics file in the workspace directory to generate the read ids for filtering.

Parsing the sequencing_summary.txt is convenient and fast. I have an old 2D example here in my test data folder, showing that it contains qscores of template, complement, and 2d. If that would be the same for 1D then filtering out 1D2 reads which we don't want to keep because they have lower quality than parents would be a breeze.

You can exclude the 1D² reads by generate a list of the ids from the 1D² reads, then filtering those reads out from the 1D fastq files. Here's one way to do that using some accessory scripts that I have written:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | awk '{print $2}' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

In order to deal with the concatenated read IDs, it's necessary to switch to a language like Perl that supports the use of backreferences for the splitting. The following perl command (replacing the first line in the above procedure) inserts a line break after the 36th character in the 2nd value, then prints that value:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt

So a full sequence of operations would look something like this:
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq
$ cat 1d2_reads.fq filtered_1d_reads.fq > unique_reads.fq


text:
QMean standalone evaluation
I want to use a script to evaluate a protein structure with QMEAN scoring.
I can't understand why the script suggested by SWISS-MODEL (access to QMEAN by json object) runs only if I put the commands on the shell. If I build a module, the result of the second request (request.get), return a QUEUEING status, not COMPLETE status.
Here is the code that I have been using:
import json
import requests

def qmean_page(structure):
""" return the Qmean4 score from SWISSMODEL QMEAN tool"""

    url = "https://swissmodel.expasy.org/qmean/submit/"

    # Request the resource
    response = requests.post(url=url, data={"email": "XXX@XXX.XX"} ,files={"structure": open(structure,'r'),},timeout=3,verify=True)

    json.dumps(response.json(), indent=4, sort_keys=True)
    download_json = str(response.json()['download_url'])

    return download_json

def qmean_score(download_url):
# Request the resource
    r = requests.get(download_url, verify=True)

    scores = json.dumps(r.json(), indent=4, sort_keys=True)
    print scores
    jdata = json.loads(scores)
    zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
    return zscore

s = 'P02489.B99990005.pdb'
url_d = qmean_page(s)
print url_d
score  = qmean_score(url_d)
print score

CASE A: If I comment the last 5 line, I can import qmean.py (the function file) in python terminal.
>>> import qmean
>>> s = 'P02489.B99990005.pdb'
>>> url_d = qmean.qmean_page(s)
>>> print url_d
https://swissmodel.expasy.org/qmean/project/TgEjha.json
>>> score  = qmean.qmean_score(url_d)
{
"input_data": {
    "info": {
        "model_001": {
            "chain_mapping": {
                " ": "A"
            }
        }
    },
    "meta": {
        "created": "2018-01-18T10:04:29.978",
        "email": "XXX@XXX.XX",
        "project_name": "Project",
        "seqres_uploaded": false
    },
    "method": "QMEANDisCo",
    "models": [
        {
            "error_status": null,
            "global_scores": {
                "acc_agreement": {
                    "norm": 0.5606936416184971,
                    "zscore": -1.0459143647477258
                },
                "all_atom": {
                    "norm": -0.012626634534455788,
                    "zscore": -2.1825427852787604
                },
                "cbeta": {
                    "norm": -0.004260555683274807,
                    "zscore": -1.9157144117276892
                },
                "qmean4": {
                    "norm": 0.6275190794418422,
                    "zscore": -3.5006904416158258
                },
                "qmean6": {
                    "norm": 0.6339071011931049,
                    "zscore": -3.0507104571131918
                },
                "solvation": {
                    "norm": -0.609099435253941,
                    "zscore": -1.2499924754236265
                },
                "ss_agreement": {
                    "norm": 0.6706229586994028,
                    "zscore": 0.05371252044245716
                },
                "torsion": {
                    "norm": 0.0033300496230023145,
                    "zscore": -2.7188811736675764
                }
            },
            "model_pdb": "https://swissmodel.expasy.org/qmean/project/TgEjha/model_001.pdb",
            "modelid": "model_001",
            "name": "P02489.B99990005.pdb",
            "seqres": [
                {
                    "atomseq": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS",
                    "chain_name": " ",
                    "name": "seq_chain_0",
                    "sequence": "MDVTIQHPWFKRTLGPFYPSRLFDQFFGEGLFEYDLLPFLSSTISPYYRQSLFRTVLDSGISEVRSDRDKFVIFLDVKHFSPEDLTVKVQDDFVEIHGKHNERQDDHGYISREFHRRYRLPSNVDQSALSCSLSADGMLTFCGPKIQTGLDATHAERAIPVSREEKPTSAPSS"
                }
            ]
        }
    ],
    "options": {
        "qmeanbrane": false,
        "qmeandisco": true
    },
    "project_type": "default",
    "sequences": []
},
"status": "COMPLETED"
}
>>>

CASE B: If I uncomment the last 5 lines and run 'python qmean.py', it doesn't work.
  andrea$ python qmean.py
  https://swissmodel.expasy.org/qmean/project/zrcLCK.json
  {
  "status": "QUEUEING"
  }
  Traceback (most recent call last):
  File "qmean.py", line 31, in <module>
  score  = qmean_score(url_d)
  File "qmean.py", line 24, in qmean_score
  zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
  KeyError: 'input_data'

If I print 'r' variable of 'qmean_score' function, It returns:
  <Response [200]>
  {
  "status": "QUEUEING"
  }

You need to use if __main__: if you want to use that file as a module and in command line. But what is your question? You seem to have solved your problems.

No, it is not solved because I could have a result like case A also when I call my function as case B

I resolve the problem by using a while loop.

maybe you could post your answer and accept it, for people trying to figure out what was the issue, or delete your question if you think that it won't be useful for other people.

So this is connected to how QMEAN is implemented as a server. It also could be true to many bioinformatics resources so I think this is important.
Posting to QMEAN "submit" triggers the job that is processed on our HPC cluster. That takes time and in the meantime JSON response just tells you that the job is QUEUEING/RUNNING or when it is ready it tells you COMPLETED. The answer suggested by you is the correct assuming you wish to wait for the job completion (which might take quite a lot).
So just check in while loop when the status is COMPLETED and fetch the result eg.
import time
import json
import requests
def qmean_score(download_url):
    # Request the resource
    zscore = None
    while True:
        r = requests.get(download_url, verify=True)

        scores = json.dumps(r.json(), indent=4, sort_keys=True)
        jdata = json.loads(scores)
        if jdata["status"] == "COMPLETED":
            zscore = jdata['input_data']['models'][0]['global_scores']['qmean4']['zscore']
            break
        else:
            time.sleep(1000)
    return zscore

But this doesn’t differ between interactive and batch usage. OP just got lucky when running the code from an interactive session.

It does in a sense that one can submit many targets, save them in eg. file and after a while, one can return to them and retrieve all the date with JSON. It is just a convenience API for people that would like to do it in the command line and not doing this by hand. It is also used in CASP and CAMEO.


text:
Looper: Duplicate sample name error
I am using Looper. While running looper run project_config.yaml I am getting duplicate sample name error. I have two fastq files for each sample, so in the annotation.csv file I have:
sample_name,library,file
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r1.fq.gz
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r2.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r1.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r2.fq.gz

My pipeline_interface.yaml file looks the following way:
protocol_mapping:
  RNA-seq: sandro_rna_seq

pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-L": "IU"
      "-T": "gencode_mouse_m13"
      "-E": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/jyvo_experiment-metadata.yaml"
      "-C": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/salmon_rna_seq.yaml"
      "-X": "True"
      "-P": "6"
      "-I": sample_name
    resources:
      default:
        file_size: "0"
        cores: "20"
        mem: "120000"

Why am I getting this error? How could I fix it?

As explained in the pep docs, the sample_name column should be considered a unique identifier for a sample. The PEP requires that One row corresponds to one sample You have multiple rows with the same name; this is not allowed.
You simply need to recast your annotation sheet such that it has file1 and file2 columns, but only 2 rows.


text:
How to retrieve SNPs data of different humans?
How can I retrieve only SNPs from the dbSNP's FTP site. The files are gzipped and organized by human chromosome (i.e, chr[1-22,X,Y]), are these for one human? I need only SNPs data for various humans.

Welcome to Bioinformatics. Could you please explain to which files are you talking (provide a link or explain the route to them)? I went to the ftp site and I couldn't find the folders. Related: have you look at this alert?

The files present at this link ftp://ftp.ncbi.nlm.nih.gov/snp/.redesign/latest_release/JSON

The readme says these files are for build 150  and the database is already by build 151

If you mean that you want snps from individuals, instead of all together, you can find them in 1000 genomes. Here different individuals from different populations are sequenced and variants are called, including snps.
Maybe you should add to your question what your goal is, what would you like to find?

Actually I am an engineer and new to bioinformatics. I am currently working on how to achieve data privacy on genomic data, for which I need human SNPs data. Like I need SNPs data for at least 100 humans. Where can I find only human SNPs data?

Did you take a look at 1000G? You can find snps of a 1000 humans!

Yes, this is the link in 1000G where I am looking ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/data/. I am not able to see any SNP sequences here. Can you please provide me the exact link if possible.

You are looking at alignments now, I guess you want the variants called already? Look for the VCF files, like here.

Thanks for the response!! I am not fully aware about the technical terms in bio. Can you please provide some information about the alignments and the variants?  How to distinguish between them?

Also are the human SNPs of same length?

Maybe you should start to read something about SNPs first, if you would know what the abbreviation SNP stands for, you would know they all have the same length of 1 nucleotide. Please read something about SNPs and next-generation sequencing first, this forum is not meant to repeat the basic information that is for example available on wiki or other websites.


text:
Looper: missing argument error
While trying to run looper I am getting the following error:

Error (missing attribute): 'sandro_rna_seq' requires sample attribute '/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-metadata.yaml' for option/argument '-E'
  Not submitted: Required attribute(s) missing for pipeline arguments string

In my pipeline_interface.yaml file I have the following:
arguments:
  "-L": "IU"
  "-T": "gencode_mouse_m13"
  "-E": "/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-
metadata.yaml"

The file "/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-metadata.yaml" is definitely present but it still gives an error...
Any suggestions would be greatly appreciated.

You are using the argument values incorrectly. Argument values are sample attributes, not full character paths. Hence, "IU" is a sample attribute, "gencode_mouse_m13" is a sample attribute, so why would you think "/scratch..." would be treated any differently?
The pipeline interface maps sample attributes to command-line arguments. It is modulated by sample (different values per sample). What you're trying to do is pass a constant string to an argument.
For this, there is a pipeline_args section of your project config. Please read the docs on project config.


text:
Splitting characters of a column in repeated lines and reducing them to one line
I have a file with sample data like this:
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042017    2   5008    T:5006  C:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     ß@  5
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042017    2   5008    T:5006  C:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     ßß  2
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042030    2   5008    C:5007  A:1     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     @@  1
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042030    2   5008    C:5007  A:1     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     @ß  3
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042029    2   5008    G:5006  T:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     @@  1
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042029    2   5008    G:5006  T:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     @ß  3
1   94853396    94853416    GCCGCCACTAGATGGTGCTA    94853398    2   5008    T:5007  C:1     94969254    94969274    CAGTCCTCTAGAGGGAGCCC    94969266    2   5008    G:5006  A:2     ß@  5
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584300   2   5008    C:5006  T:2     103841436   103841456   CATGCCACAAGAGGGCATCA    103841456   2   5008    A:5006  G:2     @ß  3
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584300   2   5008    C:5006  T:2     103841436   103841456   CATGCCACAAGAGGGCATCA    103841445   2   5008    C:5006  T:2     @@  1
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584284   2   5008    A:5007  G:1     103841436   103841456   CATGCCACAAGAGGGCATCA    103841456   2   5008    A:5006  G:2     ßß  2
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584284   2   5008    A:5007  G:1     103841436   103841456   CATGCCACAAGAGGGCATCA    103841445   2   5008    C:5006  T:2     ß@  5

As seen, some lines have repeating values for $2 and $3, line 1-6 and line 8-11 in this case.
$20 of my file contains codes of two characters representing each line, in the format "@@","ßß", "@ß","@N","ß@","ßN","N@","Nß" or "NN". Each line can have one of these 9 possible codes.
I want to write a script, which splits the characters of $20 (@,ß or N) of the lines which are repeated (1-6 & 8-11 here), and reduces them into one line in the following manner:
Count the number of @, ß or N in the first and second positions respectively of the repeated lines. For the first and second characters of each repeated line, if:

@ occurs most, print @
ß occurs most, print ß
N occurs most, print N
@ occurs equal to ß, print N
@ occurs equal to N, print @
ß occurs equal to N, print ß

In this sample dataset with repeated lines, the output file generated should be something like this, only columns $1, $2, $3, $4, $11,$12,$13, and a new column with the new "@","ß","N" codes after said counting is reported:
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90184334    90184354    GAGGCCACAAGAGGGCACAA    @N
1   94853396    94853416    GCCGCCACTAGATGGTGCTA    94969254    94969274    CAGTCCTCTAGAGGGAGCCC    ß@
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103841436   103841456   CATGCCACAAGAGGGCATCA    NN

In the first repeated case(1-6), for the first character @ occurs four times and ß occurs twice and for the second character @ and ß occur equal number of times, so the final code is @N. In the second repeated case(8-11), @ and ß for both character positions occur equally, so the final code is NN.

OK. This is a bit convoluted, but it works. There might be simpler solutions...
First we'll write a function that takes a data frame and a character position (ie; first or second position in the string), groups the data frame by the second and third columns, counts occurrences of each character in the specified position of the 18th column, and returns a vector according to your ranking scheme (this doesn't include counts of N since it doesn't occur in your example data, but it's a simple extension of the case_when line to add it):
library(dplyr)
library(tidyr)
library(stringr)

summarise_characters <- function(df, pos) {
  df %>% mutate(character=str_sub(X18, pos, pos)) %>%
  group_by(X2, X3, character) %>%
  summarise(n=n()) %>%
  ungroup() %>%
  spread(character, n, fill=0) %>%
  mutate(pos1=case_when(`@` > ß ~ '@',
                        `@` == ß ~ 'N',
                        `@` < ß ~ 'ß')) %>%
  `$`(pos1)
}

Next we remove duplicate lines from the data frame, then create a new column, pasting together the output of our function called on the first and second character:
df %>% group_by(X2, X3) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(counts=paste0(summarise_characters(df, 1), summarise_characters(df, 2))) %>%
  select(X1, X2, X3, X4, X11, X12, X13, counts)

# A tibble: 3 x 8
     X1        X2        X3                   X4       X11                  X12       X13 counts
  <int>     <int>     <int>                <chr>     <int>                <chr>     <int>  <chr>
1     1  90042011  90042031 AGTGCCACCAGGTGGGGCCG  90184354 GAGGCCACAAGAGGGCACAA  90184337     @N
2     1  94853396  94853416 GCCGCCACTAGATGGTGCTA  94969274 CAGTCCTCTAGAGGGAGCCC  94969266     ß@
3     2 103584283 103584303 TGGGCCACTAGGAGGCACTG 103841456 CATGCCACAAGAGGGCATCA 103841456     NN

This looks good. One question, what is the position input in the first function? Because I want to do this to all such repeating lines in my dataset

position in the character string in column 18. As I understand it you want to use the same logic on both the first and second position in this field

It gives me this error:  Error in mutate_impl(.data, dots) :
C stack usage  7970884 is too close to the limit

Interesting. What are the names of the columns in your data frame?

It's V1,V2,V3. But I replaced the Xs in your code with Vs already. And changed 18 to 20 since it was actually $20 as mentioned in my question.

Also I noticed that X11, X12 and X13 are in not the same order as X2, X3 and X4 as how it would be ideal. Is this causing the error?

It's column 18 in your example data. I'm not going to be able to figure out the source of the error if you're working with different data than the example

Let us continue this discussion in chat.


text:
How to use DCA to deduce the contacts between 2 proteins ensembles?
I would like to use DCA (direct coupling analysis) to analyze the interface between 2 protein ensembles (how the 2 ensembles interacts). Do you have access to a tutorial on this matter? A source code?
I found this Julia library and this library by stanford, but I am not sure if it's the right way to go and how to use these 2 for a simple test case of binding between a $\text{Ca}^{2+}$/calmodulin and a CaMKII.


text:
knnImputation for missing data
I have a dataframe containing missing values. I want to impute them with knnImputation as shown below.
Now my question is: are data imputed row wise or column wise and is there a way to give this option in knnImputation?
A <-c(1,2,3,4,5, NA, NA)
B <-c(2,3,NA, 2,4, 7, 6)
C <-c(2,3,4, NA, 5, 1,7)
D<-21:27
data<-data.frame(A,B, C, D)
new_data<-knnImputation(data, k = 3, scale = T, meth = "weighAvg", distData = NULL)

According to the code it seems like by row scaling by column but I am not sure. But you can change the orientation of your matrix by transposing it with t . Another option is to test it with a symmetric (or quasi-symmetric) matrix

The data are imputed using values from the same column
Here's a simplified example to show what's going on:
A <-c(1,2,3,4,5)
B <-c(2,3,NA, 2,4)
C <- c(2,3,4, 6, 1)
data<-data.frame(A,B,C)
data

#>  A  B C
#>1 1  2 2
#>2 2  3 3
#>3 3 NA 4
#>4 4  2 6
#>5 5  4 1

# I'm using scale = F to make things simpler
knnImputation(data, k = 3, scale = F, meth = "weighAvg", distData = NULL)

#>  A        B C
#>1 1 2.000000 2
#>2 2 3.000000 3
#>3 3 2.594272 4
#>4 4 2.000000 6
#>5 5 4.000000 1

# calculate euclidean distances, dropping column of interest
dist(data[,-2])

#>1        2        3        4
#>2 1.414214
#>3 2.828427 1.414214
#>4 5.000000 3.605551 2.236068
#>5 4.123106 3.605551 3.605551 5.099020

# use k rows with smallest distance to row of interest to calculate weighted average
# using exp(-dist) as weights
weight_sum <- weight_sum = exp(-2.828427) + exp(-1.414214) + exp(-2.236068)
2*exp(-2.828427)/weight_sum + 3*exp(-1.414214)/weight_sum + 2*exp(-2.236068)/weight_sum

#>2.594271

If you want to impute using values from the same row, all you have to do is transpose the matrix:
t(knnImputation(t(data), k = 2, scale = F, meth = "weighAvg", distData = NULL))

#>     A        B C
#>[1,] 1 2.000000 2
#>[2,] 2 3.000000 3
#>[3,] 3 3.086729 4
#>[4,] 4 2.000000 6
#>[5,] 5 4.000000 1

(I had to lower k because there are only 2 extra columns)


text:
Are mitochondrial genes to exclude in scRNA-seq such as ribosomal genes?
In this answer, it is stated that ribosomal genes should be excluded prior to normalization in scRNA-seq as contaminants.
Do mitochondrial genes have to be excluded as well? I plotted the top 50 expressed genes for a specific dataset and they tend to appear often (for example MT-ATP6). My assumption is that, given that they work for mitochondrial function and may be highly expressed, they can dilute the signal of gene differential across cell types but expressed at lower levels. Is this biologically sound?
Additionally, in this course, mitochondrial genes are used to filter cells when they contribute above a certain ratio of total RNA of single cells. However, I could not find an explanation for this filtering. What is the underlying process this procedure is used to account for?

It's worth distinguishing between the ribosomal RNA gene (Rn45s) and the many genes that code for ribosomal protein (which start with RPS, RPL, MRPS, or MRPL). mRNA's for ribosomal protein subunits are an intermediate stage used for protein production. You might catch a few lucky ribosomal RNA (Rn45s) copies that were transcribed recently and still haven't been processed, but depending on your RNA-seq protocol, there's a chance they are acting as a functional part of the ribosome rather than an intermediate stage.

According to Ilicic et al. (2016), on upregulation of mtRNA in broken cells:

There is an extensive literature on the relationship between mtDNA, mitochondrially localized proteins, and cell death [34, 35]. However, upregulation of RNA levels of mtDNA in broken cells suggests losses in cytoplasmic content. In a situation where cell membrane is broken, cytoplasmic RNA will be lost, but RNAs enclosed in the mitochondria will be retained [...]

According to this, I think that mtRNA levels relative to endogenous RNA can be used as a control for low quality (broken) cells, in a similar way to ERCC spike-ins.

Ilicic T, Kim JK, Kolodziejczyk AA, Bagger FO, McCarthy DJ, Marioni JC, Teichmann SA. Classification of low quality cells from single-cell RNA-seq data. Genome biology. 2016 Feb 17;17(1):29.

This is correct, though I'd add that while the abundance of mtRNA is indicative of broken cells, analysts should ensure that filtering done based on "percent mtRNA" should be guided by the dataset at hand, not a hard threshold.  I have seen several analysts exclude all cells with more than 5% mtRNA just because that's what is done in a Seurat tutorial.   In practice, the threshold is (a) arbitrary, (b) can vary, and (c) should be informed by visual inspection.

Mitochondrial genes do not always have to be removed, though you may need to do so for a variety of technical reasons.
Mitochondrial reads are innately different from rRNA due to rRNA typically being excluded during normal library preparation (so any of it getting carried over is essentially a contaminant). Some scRNAseq programs (e.g., RaceID) can use sampling at various stages, which can break if you have a lot of your signal going to mitochondrial transcripts. Then it's not that they should be excluded because they're from mitochondria, but rather that the results are heavily affected by the presence of VERY highly expressed genes and your results are suspect if you don't remove these RNAs. When your methods aren't affected by that (e.g., you use more robust scaling methods) or you don't have super high expression of mtRNA then that's not explicitly needed.
Further, sometimes you very much want to keep mtRNAs, since you might be interested in metabolic changes. Since I work in an institute partially dedicated to immunology, we often look at things like metabolic changes involved in immune activation and removing mtRNAs would hinder that. But if you're a priori not interested in metabolism then excluding genes you won't care about isn't unreasonable.


text:
Looper: can't specify the individual path for one pipeline run
I have multiple samples that I need to run one pipeline on. I am using Looper: https://github.com/pepkit/looper
I have a real trouble specifying a separate location for each pipeline output files. The first time when I did not care about it thinking that Looper would just automatically do that, it lead to multiple files being overwritten. Then, trying to avoid it I tried (and it did not work):

Adding {sample_name} to the output_dir of the metadata section of project_config.yaml file:

 metadata:
  sample_annotation: sandro_annotation.csv
  output_dir: /scratch/nv4e/sandro_looper_results/{sample_name}
  pipelines_dir: [../pipeline/pipeline_interface.yaml]

Creating a derived_column in the .csv annotation file:

sample_name,library,sample_name_1_read,sample_name_2_read,output_folder
abeta_24h_1,RNA-seq,read1,read2,out
abeta_24h_2,RNA-seq,read1,read2,out
abeta_24h_3,RNA-seq,read1,read2,out

pipeline_interface.yaml:
...
pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-S": sample_name
      "--input": sample_name_1_read
      "-O": output_folder
...

and project.config.yaml:
...
derived_columns: [sample_name_1_read,sample_name_2_read,output_folder]
data_sources:
  read1: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r1.fq.gz
  read2: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r2.fq.gz
  out: /scratch/nv4e/sandro_looper_results/results_pipeline/{sample_name}
...

In this particular case the parameter -O was formed correctly which I was able to see in sandro_looper_results/submissions/ .sub file, however, it was supplied the second time, apparently from the output_dir section of the project_config.yaml and was overriding the first -O command. So, obviously I tried to delete output_dir line but it is a required setting for Looper, so it can not run without it.

I have also tried supplying -O command in project_config.yaml in the pipeline_args: section with no success.

I really do not know what I could do at this point. Any suggestions would be greatly appreciated,

Looper does not generate any of these per-sample outputs.
Looper does not have any control over where your pipeline puts its outputs.
The problem is in your pipeline, not in looper. Assuming you're using a pypiper pipeline, I suggest reading the docs; for example, the "your first pipeline" tutorial: http://pypiper.readthedocs.io/en/latest/tutorials-first.html
For real-life example, see this pipeline

outfolder = os.path.abspath(os.path.join(args.output_parent, args.sample_name))
pm = pypiper.PipelineManager(name="RRBS", outfolder=outfolder, args=args, version=__version__)


text:
population admixture
I am new to Admixture and I would like to use fastNGSadmix to look for admixture in one individual (modern). Admixture is supposed to look for proportions of ancestral populations in a sample (either a population or an individual sample). for that I will need to prepare my reference population set to compare with my actual sample which is an individual form region X for instance.
Now I am thinking whether I have to add the X population to my reference set or not? If yes, wouldn't it introduce bias? Like my individual will always show the X population as the first admixing population?

Hello, welcome to Bioinformatics.SE. I have never used Admixture personally, but as far as I know it's supposed to quantify contribution of different populations to genomes (snps) of samples. So it depends on your question. You could maybe add bit more details. We like details. Details make questions answerable.

Hi, well I guess you are right about more detail. Admixture is supposed to look for proportions of ancestral populations in a sample (either a population or an individual sample). for that I will need to prepare my reference population set to compare with my actual sample which is an individual form region X for instance. Now I am thinking whether I have t add the X population to my reference set or not? if yes, wouldn't it introduce bias? like my individual will always show the X population as the first admixing population? I hope it is clear enough :)

Disclaimer, I am not an expert on admixture, answer is provided without any warranty.
It depends what you want to figure out. Why do you want to categorize the ancestry of the individual? If you analyze this origin of population X you should make a reference without population X and analyze all the individuals from the population with all the surrounding populations in the reference. However if you are interested in origin of a single individual, than the population X got to be in reference as well.
If your result will be that 90% of SNPs of your individual correspond to your population X, then it means that's a normal X fellow. Why it should be a problem?
-- edit --
I just found this preprint that highlights potential pitfalls in interpretation of STRUCTURE/ADMIXTURE results.

Thank you. well, my goal is to determine the proportions of ancestral populations and my problem with the 90% match to the source population is that I am not really sure if this 90% means that the source population in ancestral to the individual. does it? It is basically about similarities. If I compare sth with itself I will always get the best match. no?


text:
Mapping bisulfite reads to a reduced size genome
Is it worth mapping bisulfite reads (WGBS) to a reduced size genome?
I'm interested only in modifications in CpG context, thus instead of mapping to a whole genome (human genome) I would:

Extract DNA sequence around each CpG (-1000/+1000 bp)
Map reads (using bismark) to this subset

My idea is to reduce mapping time, however afraid that there might be some mapping bias or other issues.

CG sequences are pretty common; I expect you'll hit almost all of the genome when looking at a 2kb window centred on every CG, so there wouldn't be much point in doing that for efficiency purposes.
dariober has pointed out that reads coming from elsewhere in the genome could be incorrectly mapped to the reduced genome with high mapq score since the best, genuine genomic region is not present. So even if you are interested in few small regions, you should still map to the entire genome.

Don't restrict the portion of the genome you're mapping to, you'll just increase false-positive alignments by doing so. If you're worried about mapping time, then (A) use more cores or (B) use a different mapper (BSMAP was the quickest out there last I looked, though in my opinion it's fine to wait a bit longer for better quality results (I normally suggest bwameth to people)).
As pointed out by user 719016 in the comments, the above does not apply to amplicon sequencing. With that, you'll get VASTLY better performance if align to a reduced-size genome...and without an appreciably increased false-positive rate (unless your target-enrichment for sequencing is poor).


text:
Perplexity is too large
I am trying to run {Rtsne}. On pca dataset prepared by using dudi.pca from {ade4}.
For pca calculation I selected 40 components while specifying parameters for dudi.pca:
pca = dudi.pca(counts(sce_normalized), scannf=FALSE, nf=40)

{Rtsne} gives me an error: “perplexity is too large”, how could it be fixed?
The code I am trying to run is below:
name_of_file <- "dudi_pca.rds"
path_to_file <- paste0("/scratch/nv4e/scRNA-
seq/normalization/new_datasets/",
                   name_of_file, sep="")
dudi_pca <- readRDS(path_to_file)

# run Rtsne with default parameters
rtsne_out <- Rtsne(as.matrix(dudi_pca), pca=FALSE)

What are the dimensions of dudi_pca?

The dimensions are 20 000 x 40

I was using dudi_pca incorrectly. The supplied parameter to as.matrix() should have been dudi_pca$li:
rtsne_out <- Rtsne(as.matrix(dudi_pca$li), pca=FALSE)