text:
Normalize before using adonis?
I'm using adonis, an implementation of MANOVA in R (vegan package) to understand how continuous (temperature, pH) and categorical (Site ID, Month) variables affect my 16S dataset. Applying it using my metadata 'as-is', I see that ~80% of the variation (R2) is significantly explained by one variable - Site ID.
However, after suggestions from colleagues to normalize my data (as shown here for 'min-max' normalization), the same test not only shows Site ID as the third most important variable (in explaining the most variation significantly), it only shows ~31% variation being explained from the first three variables.
Any opinions on best practices for this kind of analysis?

Hi and welcome to the site. I'm afraid cross posting is very much frowned upon on the SE network. Please either delete your question from one of the two sites, or edit it so that it is specifically addressed to the subject matter of each site. For instance, for us, we would need to know what kind of data you're dealing with ("16S dataset" isn't very specific).

Also what kind of normalization did you do?  What test are you using? t-test or some package normalization and test procedure?

@Llopis that would be a "min-max" normalization, as what's found here: https://vitalflux.com/data-science-scale-normalize-numeric-data-using-r/
adonis is documented here: http://cc.oulu.fi/~jarioksa/softhelp/vegan/html/adonis.html
It uses something directly relatable to a MANOVA.

@AndréSoares Could you edit that comment into your question? Also if you could add what kind of data are you using (like terdon asked)? Why did you choose the "min-max" normalization?


text:
How to create a pathway knockout using the JavaScript library ESCHER?
I want to create my own pathway with a custom knockout function using ESCHER. Consider the pathway knockout link.
As you can see by clicking the reaction the color changes and the text with the html/css id knockout-status adds the reactions that was clicked and the growth rate.
I am trying to work out how the knockout.js  handles the clicks and changes the colors. so far I would say that this block - as commented -  handles the clicks
// initialize event listeners
var sel = builder.selection;
sel.selectAll('.reaction,.reaction-label')
    .style('cursor', function(d) {
        if (knockable(d.bigg_id)) return 'pointer';
        else return null;
    })
    .on('click', function(d) {
        if (knockable(d.bigg_id)) {
            if (!(d.bigg_id in knockouts))
                knockouts[d.bigg_id] = true;
            model = knock_out_reaction(model, d.bigg_id);
            solve_and_display(model, knockouts);
        }
    });
// grey for reactions that cannot be knocked out
sel.selectAll('.reaction-label')
    .style('fill', function(d) {
        if (!knockable(d.bigg_id)) return '#888';
        else return null;
    });

and the '.reaction,.reaction-label' finds all css class with the reactions. What I am missing is how to create my own maps with my own click listener? I checked the documentation but it handles more the creation of the maps, but I am also want to build the js knockout functionality. Could someone maybe show on a small example with 4-5 reactions how I could create my own maps.

I'm voting to close this question because it doesn't seem to relate to a specific issue, and is asking too much out of answers. Please demonstrate what you've already tried, and highlight the points where it is not working as you expect.

source: Stack Exchange

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

sort: votes


text:
Installing Rostlab's PredictProtein
I am trying to install Rostlab's PredictProtein in a Virtual Machine as detailed here https://rostlab.org/owiki/index.php/How_to_Run_The_PredictProtein_Machine_Image_with_VirtualBox:

Download PredictProtein Machine Image to your local drive. You will need the Bio-Linux-based 'VMWARE/VMDK' format version.

Choose the Debian-based version in case you do not need a graphical
  user interface. This image is only half the size of the
  Bio-Linux-based image.

Download a database image to your local drive. We recommend the 'raw' format version.
Create New Virtual Machine: Machine -> New
Name your virtual machine as 'PredictProtein'
Choose Linux and Debian as operation system
Locate and select the VMware .vmdk file as your boot disk
Give at least 2GB of RAM to your virtual machine
Click 'Finish'
Boot your new virtual machine.
If the message "The disk drive for /mnt/local-storage is not ready yet or not present" appears, press 's'.

The problem comes at this step:

Open a terminal and enter:
sudo mkfs.ext4 -F -m 0 -L local-storage /dev/sdb
sudo mount -a
sudo mkdir /mnt/local-storage/rostlab-data && cd /mnt/local-storage/rostlab-data
sudo mount -a
sudo mount /mnt/rostlab-data
sudo tar -xvJf /mnt/rostlab-data/rostlab-data.txz
sudo umount /mnt/rostlab-data && cd

When I run sudo mount /mnt/rostlab-data, I get:
mount: special device LABEL=rostlab-data does not exist

How should I proceed?
Additional notes -
The output for
cat /etc/fstab

is
# /etc/fstab: static file system information.
#
# Use 'blkid -o value -s UUID' to print the universally unique identifier
# for a device; this may be used with UUID= as a more robust way to name
# devices that works even if disks are added and removed. See fstab(5).
#
# <file system> <mount point>   <type>  <options>       <dump>  <pass>
proc            /proc           proc    nodev,noexec,nosuid 0   0
# / was on /dev/sda1 during installation
UUID=9a9418ff-48f9-421a-b25a-82161839351d / ext4 errors=remount-ro 0 1
# swap was on /dev/sda5 during installation
UUID=a96757f4-8cc7-410f-95ee-b5ae642eb9ed none swap sw  0       0

LABEL=local-storage /mnt/local-storage auto defaults,nofail 0   2
LABEL=rostlab-data /mnt/rostlab-data auto defaults,noauto,ro 0  0

/mnt/local-storage/rostlab-data /usr/share/rostlab-data none bind,nofail 0 0

The output for
blkid

is
/dev/ramzswap0: TYPE="swap"
/dev/sda1: UUID="9a9418ff-48f9-421a-b25a-82161839351d" TYPE="ext4"
/dev/sda5: UUID="a96757f4-8cc7-410f-95ee-b5ae642eb9ed" TYPE="swap"
/dev/sdb: LABEL="local-storage" UUID="df841028-8dd1-4cd2-b1bb-6ae1f6887510" TYPE="ext4"

Sounds like they've misconfigured the image. Please [edit] your question and post the output of running cat /etc/fstab and blkid from the virtual machine.

It looks like the image you downloaded is missing the partition with the label rostlab-data. Your blkid output shows that you have the local-storage partition, but not the rostlab-data which is automatically added to your /etc/fstab.
I am afraid you'll have to contact the developers and let them know. Sounds like a bug.

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


text:
Stranded vs. unstranded library preparation protocols in RNAseq
I've been reading this paper lately:
Sailfish Enables Alignment-Free Isoform Quantification from RNA-Seq Reads Using Lightweight Algorithms
I don't really understand the second paragraph under Quantification in the Online Methods section. I quote the author:

In an unstranded protocol, sequencing reads, and hence the k-mers they contain, may originate from transcripts in either the forward or reverse direction.

As far as I understand we are sequencing cDNA here and it is a single stranded molecule so all the reads will have the same direction (unlike the DNA where the reverse complement strand is also there), so I don't see the author's point. Can someone clarify please? Thanks.

Short answer: cDNA is double-stranded.
I'll try to explain the concepts surrounding this with a nanopore read, where the cDNA component of the read (from the mouse genome) is indicated in a light grey, and the adapter components are either coloured or dark grey:

The cDNA library from this RNA was prepared by adding a polyT adapter (binding to the polyA sequence of the RNA), and carrying out reverse transcription to the start of the transcript. A strand-switching process was then used to prime and extend along the opposite strand of the DNA, creating a double-stranded full-length cDNA template. This template was then amplified using PCR to increase cDNA concentration and add in barcoding sequences.
You may notice that this particular read has a polyT head, rather than a polyA tail, attached to the VNP primer in the forward orientation. I can tell by this feature that this read comes from the reverse-complement strand of a DNA transcript, and so I know that the actual RNA transcript sequence is the reverse-complement of this read. Both the forward and reverse-complement sequences could potentially come through in the sequencer, because the adapters for sequencing are attached to both sides of the cDNA template, but the transcript direction can still be determined by looking at the orientation of the adapter sequences. When mapping this sequence to the genome, only a single strand (and its associated set of kmers) needs to be considered.
When considering what happens after mapping to the genome, there are four possible match situations:

Adapter is in reverse orientation, genome is in reverse orientation
Adapter is in forward orientation, genome is in reverse orientation
Adapter is in reverse orientation, genome is in forward orientation
Adapter is in forward orientation, genome is in forward orientation

In the case of the library preparation above, situations 2 and 3 correspond to a transcript that is encoded on the forward strand of the genome, whereas situations 1 and 4 correspond to a reverse-encoded transcript. Different library preparation methods could have a different way of encoding transcript direction. See here for more details on how these match situations can be used to produce strand-specific BAM files.
If I were to use an unstranded protocol (e.g. generating a cDNA library through other means, and attaching adapters randomly to the forward or reverse transcript strand), then the adapter direction wouldn't give me information about which strand this sequence came from. I would need to consider both strands when mapping (or more correctly, when identifying the native direction of the originating transcript), because the adapter sequence has been decoupled from the transcript direction.
On Illumina systems, this strand direction calculation is performed via the chemistry on the flow cell, and the adapters are [usually] not sequenced by the sequencer at all. Adapters from one direction bind to complementary sequences on the flow cell and are extended to produce the first read, and adapters from the other direction bind to produce the second read. Even with a stranded protocol, adapters from both directions are bound (because the cDNA template is double-stranded), but the binding order (as first or second read) can be used to determine the native direction of the originating transcript.

There are a number of library preparation methods out there. Some, like the popular dUTP-based methods, result in libraries where the orientation of a read's alignment confers information regarding the strand of the transcript from whence it came. Other protocols do more PCR cycles and don't use dUTP to allow distinguishing strands. These are then termed "unstranded" libraries and one can no longer determine the strand/direction of the originating transcript from the orientation of a read's alignment. In general, a transcript will then have reads pointing in both directions in more or less equal proportions.
You tend to not see unstranded libraries much any more. Early on they were significantly more expensive, but that's not really the case any more. There also used to be two types of stranded protocols, whereas you mostly see dUTP-based protocols these days.
BTW, you'll see these same settings in most RNAseq programs. The equivalent to an unstranded protocol setting in featureCounts is -s 0, in tophat2 it was --library-type fr-unstranded (no one should use tophat2 any more, it's just an example) and so on in other programs.

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


text:
How to get the count of each kmer past 255 using khmer
I have a Fastq file and I want to get the exact count of each possible kmer from this file.
On a previous post called How to use khmer to count k-mers? Daniel Standage proposed a custom script based on khmer methods that you can see down below :
>>> outfile = open('outfile.txt', 'w')
>>> seenkmers = set()  # Consumes a lot of memory for large input!!!
>>> for read in khmer.ReadParser('reads.fq.gz'):
...     for kmer in counts.get_kmers(read.sequence):
...         if kmer not in seenkmers:
...             print(kmer, counttable.get(kmer), file=outfile)
...             seenkmers.add(kmer)

The code works if I switch counttable.get(kmer) with counts.get(kmer) but the main problem is that when I try on a customized fastq file with 300 occurrences of one kmer, the count will still be stuck at 255 and can't go further, and for some reason I can't pass through this limit.
Is this coming from the minimum table size ?
Does anyone know if there is a parameter that can work on a custom script to disable this limit ?

255 is the default maximum size of a Counttable in khmer. You want to do the following:
import khmer
counts = khmer.Counttable(31, 1e7, 1)  # Feel free to change the parameters
counts.set_use_bigcount(True)

The last line increases the maximum value from 255 to 65535. I don't think there's any way to currently go above that without changing the khmer code (at some point you have to increase the size of the number used in the hash table).

authors: Ντίνα Βαΐτση


text:
How to QC overamplified shotgun library?
I made some NEB ultra II whole-genome shotgun libraries with around ~200ng of template going into the indexing PCR step. The PCR was run for 13 cycles and I ended up having high library concentrations between 60-90 ng/µL. Should I be worried about complexity or PCR duplicates? How would all of you QC to make sure the library is complex enough and not overamplified too badly?

FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) is a good start: it will show you observed/expected ratios for the most overrepresented sequences, and running it on a few datasets you know to be high quality should give you intuition for what looks good and what doesn't.
My labmate also produced something called preseq (http://smithlabresearch.org/software/preseq/) that uses the library complexity of an initial sequencing experiment to predict how many more unique reads you would get with additional sequencing of the same library. It could give you some intuition for whether or not it is worth sequencing your library deeper to get the information you need.


text:
Creating frequency plot from weblogo with frequency unit on the Y axis
Hello I am trying to create frequency plot for a motif with frequency unit on the Y axis from 0 -100 nt as the unit as represented in this paper.(Seiler et al., 2018)

Hi @Aishwarya welcome to the site. Did you look at other questions in the site or elsewhere about making motifs? What have you tried? Which concrete problem do you have following them?

Within weblogo, add the --units probability argument and edit the figure with Illustrator (PDF/SVG) or Inkscape (SVG) to modify the y-axis, changing fractional probability ([0,1]) to a percentage ([0,100]) by multiplying by 100. Refer to the manual for more detail.


text:
Syntenic gene browser
I found that there is a Syntenic gene browser at GEvo. Do anyone know where could I find a similar browser for a local installation?

The UCSC genome browser supports local mirrors. You could make a local mirror that has the "TransMap Alignments" cross-species gene track. From the track description page:

These tracks contain cDNA and gene alignments produced by the TransMap cross-species alignment algorithm from other vertebrate species in the UCSC Genome Browser. For closer evolutionary distances, the alignments are created using syntenically filtered LASTZ or BLASTZ alignment chains, resulting in a prediction of the orthologous genes in human. For more distant organisms, reciprocal best alignments are used.

They provide a helpful tutorial on setting up a local installation here.


text:
How to retrieve the yeast locus tag from uniprot via sparql?
Let's say I want to retrieve the yeast-specific identifier for a certain protein, in the example below it will be P00330 which I would like to link to YOL086C.
When I go to uniprot's sparql UI and type
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT *
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?db up:database <http://purl.uniprot.org/database/SGD> .

}

it returns
protein db
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446

How would I now go on and find the YOL086C? When I add the line
?db ?p ?o

to the query above, I don't see anything that would directly lead me to the identifier.
I also tried to access another database
?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

but this returns all kind of databases while I am interested only in the FungiDB. How could I select this one?
How could I link the two IDs to each other?

Ok, figured it out by looking at the scheme uniprot provides:
PREFIX up:<http://purl.uniprot.org/core/>

SELECT ?protein ?yeast_id
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein up:encodedBy ?gene .
    ?gene up:locusName ?yeast_id
}

So one starts at the accession number in the middle and then uses the path to the actual locusName.
Also the path via the database works fine, I just did not use the seeAlso correctly (protein was only bound to db but not to db2, that's why it returned all entries):
PREFIX up:<http://purl.uniprot.org/core/>
PREFIX rdfs:<http://www.w3.org/2000/01/rdf-schema#>

SELECT ?protein ?db ?db2
WHERE
{
    Bind(<http://purl.uniprot.org/uniprot/P00330> as ?protein).
    ?protein rdfs:seeAlso ?db .
    ?protein rdfs:seeAlso ?db2 .
    ?db up:database <http://purl.uniprot.org/database/SGD> .
    ?db2 up:database <http://purl.uniprot.org/database/EuPathDB>

}

This will return the desired output (one would still have to process the db2 entry, of course):
protein db  db2
http://purl.uniprot.org/uniprot/P00330  http://purl.uniprot.org/sgd/S000005446  http://purl.uniprot.org/eupathdb/FungiDB:YOL086C

I am really happy to see you found your answer! Could we use this as an example query on the sparql endpoint itself?

@Jerven: feel free to use it, sure! I just updated my answer to show that also the path via the database works fine; I just used it wrong.


text:
How can I reproduce a manual NCBI search with Biopython Entrez module?
I'm trying to make a Biopython script that reproduces a manual search on the NCBI website.
My manual search gives me the following URL:
https://www.ncbi.nlm.nih.gov/nuccore?term=ANOS1[gene%20name]%20AND%20refseq[filter]
And the "Search details" box indicates:ANOS1[gene name] AND refseq[filter].
I get 338 results that I "Send to" "File" as "GI List", save as ANOS1_orthologs_manual.txt.
How can I do the same using Biopython?

This can be done by using the "Search details" as a search term in Entrez.esearch:
import os
OPJ = os.path.join
base_dir = os.getcwd()

from Bio import Entrez
Entrez.email = "my_email_here"
query = "ANOS1[gene name] AND refseq[filter]"
handle = Entrez.esearch(term=query, db="nuccore", retmax=400)
ids = Entrez.read(handle)["IdList"]
with open(OPJ(base_dir, "ANOS1_orthologs_max400_ids.txt"), "w") as ids_file:
    ids_file.write("\n".join(ids))
    ids_file.write("\n")

Rubber duck debugging at its finest :)

@terdon I did not know this expression. I learned something new today :). Regarding off-topicness, I think the question can still be useful to others, in that it gives an almost-working example of something one may try to do in a bioinformatics context.

useful, yes, but not on topic since the actual problem was a very simple mistake and therefore the answer is unlikely to help others. You could, if you want, change it around: ask how you can reproduce the manual search in python and then post your (corrected) solution as an answer. That would be great!

@terdon I see your point. I updated the question and answer.

Perfect, thanks! I would upvote, but I already had :)


text:
Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype?
So if I attempt this question using their method I put the data of genotype and phenotype into a data frame in R then use lm() to do linear regression and summary() to get an F statistic with p value 0.54. I'm not connecting what the F statistic means in this case to the answer.
Question 2: Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype? Here, we represent each SNP site as the number of minor alleles on that locus, so 0 and 2 are for major and minor homozygous sites, respectively, and 1 is for the heterozygous sites. We also assume that minor alleles contribute to the phenotype and the effect is additive. In other words, the effect from a minor homozygous site should be twice as large as that from a heterozygous site. You may use any test methods introduced in the chapter. How about permutation tests?

GENOTYPE: 1, 2, 2, 1, 1, 0, 2, 0, 1, 0
PHENOTYPE: 0.53, 0.78, 0.81, -0.23, -0.73, 0.81, 0.27, 2.59, 1.84, 0.03

Answer: Take the linear regression with F-test as an example, the p-value without the permutation test is 0.54 and is 0.55 after 1000 permutation (this number could be different with your answer because of the randomness of the permutation test).
Link to where I'm getting this: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002828#pcbi.1002828-Balding1

What part do you not understand? What did you tried? Did you look at the help pages of the test or any other resource of statistics?

Does the edit make it more clear? Sorry, the first question was poorly written

Yes, much better

Your approach is correct - it's not the only one permitted under their conditions, but it would be the standard way given only that information.
F-test compares a "full model" (predicting the phenotype using genotype and an intercept) versus a "reduced model" (predicting the phenotype using only the intercept, i.e. the mean value of phenotype). The p-value tells whether the full model predicts better, i.e. whether adding the genotype variable significantly improves the prediction. In this single-variable case, it should be approximately equal to the p-value for the beta coefficient of the genotype, if that's more familiar.
Anyway, this isn't specific to genetics, so you can check out other online sources - e.g. https://onlinecourses.science.psu.edu/stat501/node/295 for a longer, but accessible explanation of the F-test, and https://stats.stackexchange.com/ for many answers to similar questions.


text:
assigning different characters to present and absent genes from blast output using perl or python
I have file called pangenome.txt which has only one column containing the names of the genes
head pangenome.txt
    Ha1_00044161
    Ha11_00031248
    Ha16_00021664
    Ha8_00042762
    Ha4_00025807

I have another bunch of files in a tabular blast format (330 in total: individual1.blast.txt,individual2.blast.txt, ...,individual330.blast.txt ) which are the outputs of Diamond protein aligner. Below is part of individual1.blast.txt file
head blast.individual1.txt
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021664   77.4    31  7   0   94  2   49  79  7.3e-159.3
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha8_00042762    74.2    31  8   0   94  2   1   31  9.6e-158.9
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha4_00025807    77.4    31  7   0   94  2   1   31  1.3e-058.5
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00022939   77.4    31  7   0   94  2   295 325 2.1e-057.8
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha12_00032339   71.0    31  9   0   94  2   2   32  2.8e-057.4
    HWI-ST913:300:C5W5DACXX:7:1101:1477:2147    Ha16_00021553   71.0

I want to know which genes from the pangemome.txt is present in my individuals file and which are absent.
for the example above this is my desired output:
head output.individual1.txt
    Ha1_00044161   0
    Ha11_00031248  0
    Ha16_00021664  1
    Ha8_00042762   1
    Ha4_00025807   1

My desired out put is a file with 2 columns (tabular format), the first column the list of the all genes form pangemome.txt pan, and the second column is the presence and absence of these genes in the individuals files if the gene is present present in the blast.individual1.txt, I want to assign 1 and if it is absent I want to assign 0.
Can anyone help me to do this?
Thanks a lot

I wouldn't do this in Python, myself. This is a very simple text parsing problem and the standard *nix tools will be able to do it very easily. For example with awk:
$ awk 'NR==FNR{a[$1]=0; next}
       {
        if($2 in a){
            a[$2]=1
        }
       }END{
        for(i in a){
            print i,a[i]
        }
       }' pangenome.txt blast.individual1.txt
Ha1_00044161 0
Ha4_00025807 1
Ha11_00031248 0
Ha16_00021664 1
Ha8_00042762 1

The above is just a one-liner, I expanded it for clarity. You can run it on as many files as you like:
awk 'NR==FNR{a[$1]=0; next}{if($2 in a){a[$2]=1}}END{for(i in a){print i,a[i]}}' pangenome.txt blast.individual*.txt

For many large files, the fastest approach would probably be to just grep the names (unlike the awk, this will also match if the target gene is in the 1st field, but that's probably not an issue for you):
$ grep -wof pangenome.txt blast.individual1.txt.orig
Ha16_00021664
Ha8_00042762
Ha4_00025807

That will give you the list of genes found in the target blast files. To get the 1 and 0, you can do:
$ sort <(grep -owf pangenome.txt blast.individual1.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'
Ha1_00044161 0
Ha11_00031248 0
Ha16_00021664 1
Ha4_00025807 1
Ha8_00042762 1

For multiple files, just run it on all of them:
$ sort <(grep -owf pangenome.txt blast.individual*.txt) pangenome.txt |
    uniq -c | awk '{print $2,$1-1}'

Using python it should be relatively trivial, using sets. The following snippet works with your pangenome file and one blast file; adding more blast files should be trivial, making use of the union method of python's sets.
pangenome = {x.rstrip() for x in open('pangenome.txt')}
blast = {x.rstrip().split()[1] for x in open('blast.txt')}
for x in pangenome:
    if x in blast:
        print('%s\t1' % x)
    else:
        print('%s\t0' % x)

Good approach but I would use a list comprehension [] for pangenome to preserve the original order. Also you'd be safer to use the with statement to ensure the file handles are explicitly closed. There is no need to use r.strip() when defining blast. The latest python string formatting uses .format()

Actually a generator expression for pangenome would be best and most memory efficient

Here is a way to do it using a python dictionary:
#!/usr/bin/env python3

import sys

gene_filename = sys.argv[1]
blast_filename = sys.argv[2]

# Set to 0 all genes in the gene file
genes = {}
with open(gene_filename, "r") as gene_file:
    for line in gene_file:
        genes[line.strip()] = 0

# Set to 1 those genes found in the blast file
with open(blast_filename, "r") as blast_file:
    for line in blast_file:
        gene = line.strip().split()[1]
        if gene in genes:
            genes[gene] = 1

# Print the results
for (gene, presence) in genes.items():
    print(gene, presence, sep="\t")

exit(0)

Executing it:
$ ./find_genes_in_blast.py pangenome.txt blast.individual1.txt
Ha1_00044161    0
Ha11_00031248   0
Ha16_00021664   1
Ha8_00042762    1
Ha4_00025807    1

thanks for the code. I tried it, but I get an error message:   File "./find_genes_in_blast.py", line 23
print(gene, presence, sep="\t")
^
SyntaxError: invalid syntax

If you are using python 2, you need to put from __future__ import print_function at the beginning of the code. Or you need to use a correct syntax for print in python 2, for instance print "%s\t%d" % (gene, presence).


text:
How to read .bwt file .sa file?
In BWA project in bwt.c file I found bwt_dump_bwt method and bwt_dump_sa method. I want to use them to read .bwt file and .sa file. I wrote the following program:
int main(void)
{
  char *tmp, *prefix;
    bwt_t *bwt;
    prefix = bwa_idx_infer_prefix("hg19bwaidx");

    tmp = calloc(strlen(prefix) + 5, 1);
    strcat(strcpy(tmp, prefix), ".bwt");
    bwt = bwt_restore_bwt(tmp);

  return 0;
}

I got the following error:
Undefined symbols for architecture x86_64:
  "_bwa_idx_infer_prefix", referenced from:
      _main in readbwt-3022ec.o
  "_bwt_restore_bwt", referenced from:
      _main in readbwt-3022ec.o
ld: symbol(s) not found for architecture x86_64
clang: error: linker command failed with exit code 1 (use -v to see invocation)

Any help ?
Thanks.

In the bwa source code, file bwt.{c,h} implement basic operations on .bwt and .sa files. They depend on utils.{c,h}, but you can create a fake utils.h to cut this dependency as is shown below. Then you only need three files: kvec.h, bwt.h and bwt.c to access .bwt and .sa.
The following example demonstrates how to find super-maximal exact matches (SMEMs). Note that internally bwa concatenates all database sequences into one long string. If you use bwt.{c,h} only, you get coordinates on this long string but not on each chromosome. You have to call APIs in bntseq.{c,h}, which are fairly independent, too. You can go on to implement a mini bwa-mem with this example, but you will need all the bwa-mem source code, except the main function.
Compiling and running
# compile the example
git clone https://github.com/lh3/bwa
mkdir -p tiny-smem && cd tiny-smem
cp ../bwa/{bwt.c,bwt.h,kvec.h} .              # you need at least these three files
### then create utils.h and test.c shown below ###
gcc -g -O2 -Wall bwt.c test.c -o tiny-smem    # compile tiny-smem

# create a bwa index if you don't have it. We use chr20 from GRCh38 as an example:
(cd ../bwa && make)                           # compile bwa for indexing
../bwa/bwa index -p chr20 path/to/human/chr20.fa

# run the example
./tiny-smem chr20.bwt chr20.sa TTTTCTTCTTTTGTAATATAAGCATTgGTTGTTCTAAATTTCTCTCTTGGCACTGCTTTA

The last command line above outputs:
0       0       26      1       -125734
0       14      27      4       -3657204        +16021887       +11802682       +54692983
0       15      29      1       +8344885
0       17      31      1       -39545379
0       19      37      1       +55060482
0       25      39      2       +1198672        -26409692
0       27      60      1       -125700

where column 2 and 3 give the coordinate of SMEM on the query, col4 the number of exact hits and the following columns are the list of hits in the concatenated coordinate.
File: utils.h
#ifndef FAKE_UTILS_H
#define FAKE_UTILS_H

#define xassert(cond, msg) assert(cond)
#define xopen(fn, mode) fopen((fn), (mode))
#define err_fclose(fp) fclose(fp)
#define err_fread_noeof(ptr, size, n, fp) fread((ptr), (size), (n), (fp))
#define err_fwrite(ptr, size, n, fp) fwrite((ptr), (size), (n), (fp))
#define err_fflush(fp) fflush(fp)
#define err_fseek(fp, off, whence) fseek((fp), (off), (whence))
#define err_ftell(fp) ftell(fp)

#endif // ~FAKE_UTILS_H

File: test.c
#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include <stdio.h>
#include "bwt.h"

int main(int argc, char *argv[]) {
  bwt_t *bwt;
  int j, max_intv = 100; // ignore SA intervals larger than this value
  bwtintv_v r = {0,0,0};

  if (argc < 4) { fprintf(stderr, "Usage: ./test <.bwt> <.sa> <seq>\n"); return 1; }
  if ((bwt = bwt_restore_bwt(argv[1])) == NULL) abort(); // read .bwt
  bwt_restore_sa(argv[2], bwt);                          // read .sa
  for (j = 3; j < argc; ++j) {
    int i, k, x, len = strlen(argv[j]);
    uint8_t *seq = (uint8_t*)malloc(len);
    for (i = 0; i < len; ++i) { // convert A/C/G/T/N => 0/1/2/3/4
      int c = tolower(argv[j][i]);
      if (c == 'a') seq[i] = 0;
      else if (c == 'c') seq[i] = 1;
      else if (c == 'g') seq[i] = 2;
      else if (c == 't') seq[i] = 3;
      else seq[i] = 4;
    }
    x = 0, r.n = 0;
    while (x < len) {
      x = bwt_smem1(bwt, len, seq, x, 1, &r, 0); // find SMEMs
      for (k = 0; k < r.n; ++k) { // traverse each SMEM
        bwtintv_t *p = &r.a[k]; // p->info encodes start and end on the query
        printf("%d\t%d\t%d\t%d", j-3, (int)(p->info>>32), (int)p->info, (int)p->x[2]);
        if (r.a[k].x[2] <= max_intv) { // don't retrieve coordinates for huge intervals
          int l_mem = (int)p->info - (int)(p->info>>32);
          for (i = 0; i < p->x[2]; ++i) { // traverse each hit
            uint64_t y = bwt_sa(bwt, p->x[0] + i); // get suffix array coordinate
            if (y < bwt->seq_len>>1) printf("\t+%ld", (long)y);
            else printf("\t-%ld", (long)(bwt->seq_len - y - l_mem));
          }
        }
        putchar('\n');
      }
    }
    free(seq);
  }
  free(r.a);
  bwt_destroy(bwt);
  return 0;
}

You need to include bwa.o in your call to clang, since it holds the definitions of those functions. This will cause a good number of object files, so you might as well just clang -o foo *.o foo.c -lz -lpthread -lm -lrt and get it over with. Since this was quite obviously never meant to be done, you'll have to remove at least one main() function from bwa.c.
Note that further questions regarding this will be considered off-topic and closed, since it's purely about compiling something that you want to write yourself.

Thanks. But there is no main() in bwa.c by the way I want to do something ppl did before. like this call to bwt_dump_bwt in this link: https://github.com/nsoranzo/bwa/blob/master/bwtindex.c The problem is I'm new to c programming. That why I need a help and still I can't understand how to do it.

There's a main in a different file then. You linked to a fork of the bwa repository, which apparently used to include that file.

I ran clang -o main bwa.o main.c -lz -lpthread -lm -lrt and I got: ld: library not found for -lrt and when I removed  -lrt I got many other errors. Any idea? I'm a Mac user maybe this is the problem.

You'll have to figure out the equivalent library (or install it) on OS X, or whatever system you're using.

I ran BWA successfully on my system. But I can't read .bwt with the program that I wrote.


text:
bed file merge book-end features only if same name in column 4
I have a bed file where some entries have the same name in column 4 but are book-ended in coordinates, e.g. see below:
chr1    0       10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   12000   DDX11L1::chr1:11868-14409.E1
chr1    12000   12200   DDX11L1::chr1:11868-14409.E1
chr1    12200   13000   DDX11L1::chr1:11868-14409.E1
chr1    13000   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   15800   WASH7P::chr1:14403-29570.E1
chr1    15800   16000   WASH7P::chr1:14403-29570.E1
chr1    16000   18200   WASH7P::chr1:14403-29570.E1
chr1    18200   18400   WASH7P::chr1:14403-29570.E1
chr1    18400   23000   WASH7P::chr1:14403-29570.E1
chr1    23000   23400   WASH7P::chr1:14403-29570.E1
chr1    23400   24000   WASH7P::chr1:14403-29570.E1
chr1    24000   24200   WASH7P::chr1:14403-29570.E1
chr1    24200   25200   WASH7P::chr1:14403-29570.E1
chr1    25200   25800   WASH7P::chr1:14403-29570.E1
chr1    25800   26200   WASH7P::chr1:14403-29570.E1
chr1    26200   26400   WASH7P::chr1:14403-29570.E1
chr1    26400   27000   WASH7P::chr1:14403-29570.E1
chr1    27000   27600   WASH7P::chr1:14403-29570.E1
chr1    27600   27800   WASH7P::chr1:14403-29570.E1
chr1    27800   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   30800   MIR1302-2HG::chr1:29553-31109.E1
chr1    30800   31400   MIR1302-2HG::chr1:29553-31109.E1

I would like to apply an operation that seems to me is a mix of bedtools merge -d 0 and bedtools groupby, but I can't figure out what's the combination.
For example, the 4 entries with DDX11L1::chr1:11868-14409.E1 I would like merged as they are book-ended, same with the WASH7P::chr1:14403-29570.E1 entries and the MIR1302-2HG::chr1:29553-31109.E1 entries. There are two LOC101929784::chr1:13201-13800::E1 entries, but I don't want to merge them because they are not book-ended, since they are separated by LOC101929784::chr1:13201-13800::XYZ.
Any ideas?

awk 'BEGIN{OFS="\t"}{if($1==lchrom && $4==lname && $2 == lend) {lend = $3}else{if(lchrom) {print lchrom, lstart, lend, lname;}; lchrom=$1; lstart=$2; lend=$3; lname=$4}}END{print lchrom, lstart, lend, lname}' foo.bed > new.bed

That produces:
chr1    0   10000   LOC101929784::chr1:13201-13800::E1
chr1    10000   10600   LOC101929784::chr1:13201-13800::XYZ
chr1    10600   10800   LOC101929784::chr1:13201-13800::E1
chr1    10800   13200   DDX11L1::chr1:11868-14409.E1
chr1    13200   29600   WASH7P::chr1:14403-29570.E1
chr1    29600   31400   MIR1302-2HG::chr1:29553-31109.E1

This explicitly checks for bookends and nothing else and assumes that everything is sorted.


text:
Combine multiple species trees
I have five species trees run using starBeast and Beauti. Here is a link to the paper discussing both,
http://journals.plos.org/ploscompbiol/article/file?type=supplementary&id=info:doi/10.1371/journal.pcbi.1003537.s006
 Each of these trees has the same individuals, and each locus used (10 loci per tree) has individual models set to it. All the trees used a relaxed log clock normal. I would like to be able to combine these trees into one. I have attempted to find a way to do this within Densitree and starBeast, but to no avail. When all loci are run with all trees, the program simply dies, leading me to believe that running the trees together will all parameters is computationally impossible for my hardware. Is there a way to combine these trees?

I am not familiar with trees but what are these software that you use? Could you provide a link to see what are they capable? What is a "relaxed log clock normal"? Why don't you use all 30 loci to build a single tree by combining the individual models? What do you mean that is computationally impossible: how much time do you expect to take (with how many cpus?) ?

Although I've never tried it myself, I think what you need can be achieved through the LogCombiner program from the BEAST package -- assuming the species trees files are comparable (ran with same parameters, same species names, etc.) I am also assuming that you meant 'five species tree files', each file with the posterior sample of species trees.
However, I do believe you should try to run all loci within the same analysis, or at least do some prior check to see if the distributions between data sets are comparable. LogCombiner assumes independent runs from the same data.
I would also recommend taking a look (and following) the BEAST google group (which includes starBEAST discussions).


text:
bed file with N regions of GRCh38 reference?
How do I obtain a bed file with the list of N-regions of GRCh38 reference? This is, the regions where the sequence is a stretch of Ns.

# if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk && make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed

Option -n sets the min stretch length. Just use -p as is. It is a bit complicated to explain what it is doing.

This information is stored in the 2bit file representation of the sequence, so if you happen to have a 2bit file locally (or want to download one from UCSC) and have py2bit installed (you'll need version 3.0, since I literally just added support for this):
import py2bit
tb = py2bit.open("genome.2bit")
of = open("NNNNN.bed", "w")  # Do change the input and output file names
for chrom in tb.chroms().keys():
    blocks = tb.hardMaskedBlocks(chrom)
    for block in blocks:
        of.write("{}\t{}\t{}\n".format(chrom, block[0], block[1]))
of.close()
tb.close()

In case it's helpful, you can also use this to get all soft-masked regions. Just replace hardMaskedBlocks() with softMaskedBlocks() and ensure to specify storeMasked=True when you open the file.

Use twoBitInfo:
$ twoBitInfo file.fa -nBed output.bed

For example, to get all the N-masked regions on chromosome Y (also note you can use stdout as a filename to write directly to stdout, and use of the url as the input, no need to download the 2bit file):
$ twoBitInfo http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.2bit:chrY -nBed stdout | head -3
chrY    0   10000
chrY    44821   94821
chrY    133871  222346

You can download twoBitInfo from the directory appropriate to your operating system here:
http://hgdownload.soe.ucsc.edu/admin/exe/

I might as well jump on as well. Here's a Perl script that I wrote a while ago to split a fasta sequence at stretches of Ns:
https://github.com/gringer/bioinfscripts/blob/master/fasta-nsplit.pl
I've just modified it to spit out a BED formatted file showing the N locations on standard error. Use it as follows:
./fasta-nsplit.pl tig87_withN.fa 2>out.bed > out.split.fa

Output (BED file):
tig00000087 0   60
tig00000087 900 960
tig00000087 3840    3960
tig00000087 14880   14940
tig00000087 59520   59700
tig00000087 93000   93120
tig00000087 107880  107940
tig00000087 109135  109140

Here's the full script, for completeness / fixing / comment:
#!/usr/bin/perl
use warnings;
use strict;

my $seq = "";
my $shortSeqID = "";
my $seqID = "";
my $keep = 0;
my $cumLength = 0;
while(<>){
  chomp;
  if(/^>((.+?)( .*?\s*)?)$/){
    my $newID = $1;
    my $newShortID = $2;
    if($seq){
      my $inc = 0;
      while($seq =~ s/(NNNN+)(.*)//){
 my $nStretch = $1;
        my $newSeq = $2;
        printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
 $cumLength += length($seq);
 printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength,
        $cumLength + length($nStretch));
 $cumLength += length($nStretch);
        $seq = $newSeq;
      }
      printf(">%s\n%s\n", $seqID, $seq) if ($seq);
    }
    $seq = "";
    $shortSeqID = $newShortID;
    $seqID = $newID;
    $cumLength = 0;
  } else {
    $seq .= $_;
  }
}
if($seq){
  my $inc = 0;
  while($seq =~ s/(NNNN+)(.*)//){
    my $nStretch = $1;
    my $newSeq = $2;
    printf(">%s.%s\n%s\n", $seqID, $inc++, $seq) if ($seq);
    $cumLength += length($seq);
    printf(STDERR "%s\t%d\t%d\n", $shortSeqID, $cumLength,
       $cumLength + length($nStretch));
    $cumLength += length($nStretch);
    $seq = $newSeq;
  }
  printf(">%s\n%s\n", $seqID, $seq) if ($seq);
}

Here's a way to generate it yourself from the genome sequence. First, convert the fasta file of the genome to tbl format (<seq id>\t<sequence>), then use perl to find the start and end positions of all stretches of consecutive of N or n.
FastaToTbl hg38.fa |
    perl -F"\t" -ane 'while(/N+/ig){
                            for(0..$#-){
                                print "$F[0]\t$-[$_]\t$+[$_]\n"
                            }
                      }' > hg38.n.bed

Explanation

FastaToTbl : this is a very simple script that converts fasta to tbl. Just save the lines below somewhere in your $PATH (e.g. ~/bin) as FastaToTbl and make it executable.
#!/usr/bin/awk -f
{
    if (substr($1,1,1)==">")
            if (NR>1)
                printf "\n%s\t", substr($0,2,length($0)-1)
            else
                printf "%s\t", substr($0,2,length($0)-1)
        else
            printf "%s", $0
}
END{printf "\n"}

The Perl magic.

The -a makes perl behave like awk, splitting each input line on the value given by -F (a tab in this case) and saving the result in the aray @F. So, $F[0] will be the sequence id.
while(/N+/ig){ : this will match all cases of consecutive N or n (the i flag makes it case-insensitive). Perl will store the start positions of all matches in the special array @- and the end positions in @+.
for(0..$#-){ : iterate over all numbers from 0 to the final index ($#-) in the array @-.
print "$F[0]\t$-[$_]\t$+[$_]\n" : print minimal bed format data. The name of the current sequence ($F[0]), the start position of the current match ($-[$_]) and the corresponding end position ($+[$_]).

I just ran the above on my system and it generated this.

@user172818 code can't count N at fa begining
@terdon code get wrong bed position
i change @terdon code like below; and it works for me

awk '{
    if (substr($1,1,1)==">")
        if (NR>1)
            printf "\n%s,", substr($0,2,length($0)-1)
        else
            printf "%s,", substr($0,2,length($0)-1)
    else
        printf "%s", $0
}
END{printf "\n"}' test.fa | \
perl -F"," -ane 'while($F[1] =~ /N+/ig){
                        for(0..$#-){
                            print "$F[0]\t$-[$_]\t$+[$_]\n";
                        }
                  }'

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


text:
I have two files: mys.nuc and ms.prep and need to use bioseq to extract the ORF sequences and save as “mys.nuc2”
Below are the two files I am working with:
$ cat mys.pep
>00001 16 324  len=309
LDRFYTKKVIMYKQQYFISGKVQGVGFRFFTEQIANNMKLKGFVKNLNDGRVEIVAFFNT
KEQMKKFEKLLNGNKYSNIKNIEKIVLDENYPFQFNDFKIYY*
>00002 751 308  len=398
YKLYNILLIKYFLIFPKYIDNA*I*KNKNYFKTFLLDRFYTKKVIMYKQQYFISGKVQGV
GFRFFTEQIANNMKLKGFVKNLNDGRVEIVAFFNTKEQMKKFEKLLNGNKYSNIKNIEKI

$cat mys.nuc
>00001  16 324  len=309
TTGGATAGATTTTATACAAAGAAGGTAATAATGTATAAACAACAATATTTTATTTCTGGC
AAGGTGCAAGGTGTTGGTTTTAGATTTTTCACAGAGCAAATAGCAAATAATATGAAACTA
AAAGGATTTGTAAAAAATCTCAACGATGGAAGGGTAGAAATTGTAGCTTTCTTTAATACT
AAAGAACAAATGAAAAAATTTGAAAAATTATTAAATGGGAATAAGTATTCAAACATTAAA
AACATTGAAAAAATAGTTTTAGATGAAAATTATCCTTTTCAATTTAATGATTTTAAAATT
TATTATTAG
>00002  751 308  len=398
TACAAATTATATAATATCTTATTAATAAAATATTTCCTAATATTTCCCAAATATATTGAT
AATGCCTGAATTTAAAAAAACAAAAACTATTTTAAAACTTTTCTATTGGATAGATTTTAT
ACAAAGAAGGTAATAATGTATAAACAACAATATTTTATTTCTGGCAAGGTGCAAGGTGTT
GGTTTTAGATTTTTCACAGAGCAAATAGCAAATAATATGAAACTAAAAGGATTTGTAAAA
AATCTCAACGATGGAAGGGTAGAAATTGTAGCTTTCTTTAATACTAAAGAACAAATGAAA
AAATTTGAAAAATTATTAAATGGGAATAAGTATTCAAACATTAAAAACATTGAAAAAATA
GTTTTAGATGAAAATTATCCTTTTCAATTTAATGATTT

What I have done so far is try to run:
$ bioseq -B mys.nuc

and get:
00001
00002

Ultimately I am mean to extract the ORF sequence and save as mys.nuc2 and then translate mys.nuc2 with bioseq -t1 where that file should NOT contain any internal stop codons.
I do not understand why I get:
00001
00002

Or what I am meant to do with it, if that indeed is the correct output.

Can you run, which bioseq and whereis bioseq?  Also could you add the head bioseq to be sure? How confident are you with writing your own script to find ORF, or download and install (if you can) new tools to do it for you?


text:
Combine Fastq by writeFastq is not working properly
Currently I have 4 fastq file for each patient as WES data. For WES pipeline, I want to combine those fastq files. I am using writeFastq from ShortRead package in Bioconductor.But I think this function is not working properly.
The reason of my doubt is the the size the next step is quality check using FastQC tool. This tool is working on individual file but not on combined file and shows the error that ID line doesn't start with "@" . So I want to know what are the other effective tool for combining multiple fastq files without any loss of information.
Let me know if I am doing something wrong. Thanks
[EDIT-1] : The codes for combining multiple fastq files into one is as follows:
reference_fastq = file.path(paste(getwd(),"refs",sep = "/", collapse = " "), fsep = ".fasta")
fastqDir <- file.path(dataDir, fsep = ".fastq")
outfile_temp <- paste(dirname(dataDir),basename(fastqDir),    paste(basename(fastqDir),"combined","fastq", sep = "."),sep = "/", collapse = " ")
fls <- list.files(fastqDir, "fastq$", full=TRUE)
for (fl in fls) {
  fq <- readFastq(fl)
  writeFastq(fq,outfile_temp, mode="a")
}

First four lines of my fastq file (one of 4 indivual file) are as follows :
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

First four lines of combines fastq file are as follows:
@D00239:185:C9KELANXX:2:2201:1088:2107 1:N:0:TCCTGAGC
GTACAAGAAGCATAACATCAGCATGTGTTTTGGTGAGAACCTTAGGAAGTTTGCAATCATGACAGAAGGTGAAGG
+
/<BBBFFFBFFFFFFFFFFFBF<FFBF<B<FFFFFFFFFFB/B/<////<<<FF<B<FF<FFFFFFFFF/FFFFF

It is pretty clear that combined fastq file is also starting wiht proper identifier but why FastQC tool is saying that identifier '@' is missing. Here is the screen shot of the error message:

Thanks.

Can you please show the first and last lines of the input FASTQ files, and the first and last line of whatever writeFastq has created?

Could you paste the command you used from ShortRead to combine the files? (And the version of the package, just in case there is a bug in that version). If you have the fastq files compressed you can simply cat them to a new .tar.gz file and you'll have them combined

@Lot, please instead of images paste the text in code formatting. Also could you, answer to gringer?

@gringer: I have shown the first 4 lines of input fastq file and output fastq file. What else you want to check..let me know.

@Llopis: Actually I tried to copy the text instead the image but I am working on server and on right click after text selection, I did not find any option like copy. That's why I took screen shot. Sorry for this.

@Llopis: ShortRead package version is 1.28.0, R version 3.4.3 and bioconductor version is 3.2.

You have a 2 year's old version of Bioconductor, you might be interested in upgrading it (just in case some other errors are already solved in the updated version).

You can drastically simplify your effort because FASTQ files are simple text files; therefore, standard UNIX test file tools work:
cat refs/*.fastq > combined.fastq

And this even works with .fastq.gz files — no unzipping required.

Apart from this, your R code contains several errors. For instance, your variable names change throughout, and you’re using the fsep parameter wrongly. In fact it’s unclear to me what the code is supposed to do. At an approximation (and untested), the following should work:
fastq_dir = file.path(getwd(), "refs")
fastq_files = dir(fastq_dir, pattern = "\\.fastq$", full.names = TRUE)
outfile = file.path(dirname(fastq_dir), "combined.fastq")

for (file in fastq_files) {
    fq = readFastq(file)
    writeFastq(fq, outfile, mode = "a")
}

This is equivalent to running the shell command shown above, just in a more convoluted way (parsing and re-generating the FASTQ files).

@Konard : Thanks for pointing out my mistake. As I am beginner so my codes are not as much good. I have changes this. Can you tell me what is the proper use of fsep. I though using file path, fsep is the way to choose only specific type of file (may be I am wrong). And how can I use cat within for loop because it is giving unexpected symbol error by trying cat fq > combine.fastq Thanks again. I'll try my best to improve. keep pointing my mistakes. :)


text:
Human Cell Atlas - Data availability
A news item from 2017-10-18 on the website of the Human Cell Atlas states:

In addition, the consortium today also announced the impending release of gene expression profiles from the first one million immune cells collected under the HCA, toward an initial milestone of collecting at least 30 million cells representing several tissues and organs for the atlas’ first draft. These data, to be posted on an online repository by early November, will be freely available for researchers’ use. (https://www.humancellatlas.org/news/14)

As of now, Feb 2018, I cannot find these data anywhere.
Does anybody know where I can download the gene expression data of these 1 million immune cells or when it will be eventually published?

The preview version is now available here

There is now a post on Biostars:

We are starting to define the release process for the Human Cell
  Atlas.
The first couple of HCA releases will be primarily single-cell RNA-Seq
  and we are interested to know what you would expect our releases to
  contain
Please fill out our survey if you would like to help us define our
  releases
https://www.surveymonkey.com/r/S6Y753P

And I suppose that answers the original question: the data has not been released yet and probably won't be for a while considering they are still discussing how to make it available.
Update (from the same post):

The first production release of the DCP is happening in June.
The first data release of HCA is targetting Q4 2018.


text:
How do programs like STRUCTURE determine populations from alignments?
If I were to run BAM/SAM files through a program like STRUCTURE, or ADMIXTURE described  here how would they determine which individuals belong to what groups?
What part of the files would these programs use to determine which individuals belonged to different populations? I believe they do this via determining allele frequencies for each individual and then putting them into a population based on allele frequency, is the only criteria that they use comparison of allele frequencies between individuals?

structure works with SNP/genotype information, so won't work directly off BAM/SAM files. You would need to convert to VCF (e.g. using mpileup/vcftools), then process using something like PLINK to convert it into a format that STRUCTURE recognises.
Here's the description I have about structure from my PhD thesis. It might be a bit dated now, but I expect most of it is still applicable:
A clustering program called structure was created by Falush et al. (2003) that is used to determine population structure using a set of unlinked genotypes. The expected use case for structure is in the analysis of population structure, but it can also be useful for determining the utility of a given set of markers for categorising (or quantifying) a given trait.
The main benefit of using structure over other methods is that it is able to derive estimated ancestry coefficients from the data, rather than a discrete yes/no for each individual. This allows for the use of cutoff values for the assignment of individuals (useful for tweaking false positive and false negative rates in diagnostic tests), and a quantifiable value where a classification is not useful or does not make sense (as in the case of a continuous trait with high heritability such as height).
The program begins by assigning each individual to a given population. This can be either given in the input data as a ‘popinfo’ flag, or determined randomly by structure. After this, the program constructs a probability distribution for the data, considering the population assignments given, and uses that to estimate population allele frequencies. The program then constructs another probability distribution for these data, and uses that to estimate the population that each individual is likely to have originated from. This process is repeated many times (such a process is known as a random walk), and the underlying theory suggests that the progression of population assignments through the random walk process can be used as an approximation of the actual populations of origin for each individual.
The documentation included with structure describes the general algorithm for each iteration (of maybe hundreds of thousands of iterations) as follows:

Sample (population allele frequencies for iteration m) from the probability of (allele frequency) given (genotypes of sampled individuals) and (population of origin of individuals for iteration m − 1).
Sample (population of origin of individuals for iteration m) from the
probability of (population of origin) given (genotypes of sampled individuals) and (allele frequencies for iteration m).

In more simple terms, the steps are similar to the following statements:

these populations probably have this allele frequency distribution
therefore...
these individuals probably came from this population

The structure program is able to accommodate admixture into its model
by considering the population of origin as a probability, rather than a single
discrete value.
The file that is output from structure is a text file with individual identifiers and predicted ancestry coefficients (Q values). While the program does provide means to display these data graphically, the visual representation of data (including sorting) is better controlled using an external program. The most commonly used program is Noah Rosenberg’s distruct (Rosenberg, 2004), a program that allows colour customisation and group re-ordering. For the purpose of this thesis, a custom R script, snpchip2structure.R, has been designed which also allows selection of different sorting methods, a scatter plot for two populations, error bars, and a few additional features (see figure 1.7). More details about this script can be found in the programs directory of the supplementary CD for this thesis.


text:
NPDock - service unavailable
NPDock service http://genesilico.pl/NPDock/ is unavailable, at least that's the message that shows after you try to submit a job. It's been like this for the last few weeks. Contact email address provided on the page also doesn't work. Does NPDock work for anyone else?

I also got the same error message; the website also says "Please run max. 5 jobs at the same time. If you need to run more jobs, please consider installing NPDock (DARS-RNP) locally.", but I couldn't find any link to the actual software.

I'm voting to close this question as off-topic because it is a question that should be asked of the tool's developers and not a general bioinformatics issue.

Sorry but we can't help with this sort of thing.

Ok, now it's fixed.


text:
Why can using more reads lead to a lower quality assembly?
I am experimenting with adding additional reads to the input files I'm giving SOAPdenovo2, and there comes a point where a good contig I've been watching actually stops showing up. Does anyone have a quick answer to why that might happen?

Are you doing any sort of digital normalization or correction of the input reads, or just feeding increasing numbers of raw reads in?

They are reads that have been trimmed for quality, and have been selected from a larger read pool based on mapping the pre-assembled contigs.

Excessively high depth may introduce errors that can't be fixed. Some assemblers subsample reads to a lower depth. In general, de novo assembly is a mess. Even assembler developers often can't explain why things happen. Let alone outsiders.

The more reads you add the more errors you add into the assembly. This is because additional duplicate reads don't add nodes/edges to the de Bruijn graph, but those with errors do. By preselecting those reads aligning to your desired contigs, you may be further exacerbating this effect and further hindering SOAPdenovo2's ability to do k-mer correction (after all, the observed k-mer frequency of errors would be higher). This then results in additional incorrect paths through the graph that the assembler has to traverse.

Yes, this is commonly observed and both the previous answers (Devon Ryan and conchoecia) are spot on so I won't reiterate their explanations on why.
Try bfc (bloom-filter correction) to "fix" likely erroneous base-calls (i.e. uncharacteristically rare kmers) in the reads prior to assembly.  This will not only help your assembly, but it shall require less RAM and time.
I also recommend trying running your assembler with multiple kmer parameters, as some values will work better than others for a particular dataset.  You can simply choose the best assembly or go through a more complicated effort to merge them, depending on your circumstances.

I think the reason that no one has submitted an actual answer yet is because this question is really dependent on the size of the genome you are sequencing, its repetitive content, the level of heterozygosity, the amount of contamination, the type of data (illumina, long illumina, long reads), the library type, library peculariaties, PCR error, the heterozygosity of the contig in question, the complexity of the contig in question (high information content? repeats? AT/GC%?), and probably some other things that I haven't thought of.
Like @user172818 said, all of those things ON TOP OF all of the quirks of different assemblers make this a hard question to answer for you without a full diagnosis of your data and the genome of the thing you are sequencing/assembling.


text:
Alignment QC differences between HISAT2 and Qualimap
I used hisat2 for aligning reads to to the genome. I have an alignment summary for sample1 as follows:
89038751 reads; of these:
  89038751 (100.00%) were paired; of these:
    5641019 (6.34%) aligned concordantly 0 times
    67665552 (76.00%) aligned concordantly exactly 1 time
    15732180 (17.67%) aligned concordantly >1 times
    ----
    5641019 pairs aligned concordantly 0 times; of these:
      651415 (11.55%) aligned discordantly 1 time
    ----
    4989604 pairs aligned 0 times concordantly or discordantly; of these:
      9979208 mates make up the pairs; of these:
        5887822 (59.00%) aligned 0 times
        2656902 (26.62%) aligned exactly 1 time
        1434484 (14.37%) aligned >1 times
96.69% overall alignment rate

To check the alignment quality I used Qualimap. In the report I see the following:
Reads alignment

Number of mapped reads: 172,189,680
Total number of alignments: 244,302,519
Number of secondary alignments: 72,112,839
Number of non-unique alignments:    22,973,098
Aligned to genes:   2,081
Ambiguous alignments:   173
No feature assigned:    178,975
Not aligned:    5,887,822

Reads genomic origin

Exonic: 2,081 / 1.15%
Intronic:   6,866 / 3.79%
Intergenic: 172,109 / 95.06%
Intronic/intergenic overlapping exon:   302 / 0.17%

Transcript coverage profile

5' bias:    0.69
3' bias:    2.85
5'-3' bias: 0.39

Can anyone tell me why the numbers look different? Did I go wrong somewhere? Help me in understanding the results.

89038751 pairs times 2 is 178 million reads. That times the approximately 97% alignment rate is around 172 million reads, which is what qualimap is reporting.
A "non-unique" alignment is one that can map to multiple places (either equally well, or where all of the alignments are of similar enough quality). Such alignments have a single "primary" alignment and one or more "secondary" alignments, which is where the "secondary" metrics come from.
For the qualimap output, there are additional metrics describing how reads overlap genes (or not). Of note are the "ambiguous" alignments, which are those overlapping multiple genes. "No feature" metrics, then, refer to alignments not overlapping any gene (see the qualimap documentation for further details).
My guess is that you have such a high "intergenic" rate because you've either used the wrong GTF file or have told qualimap the wrong strand setting. For unusual organisms one expects many reads to arise from unannotated genes. This is not the case for common model organisms, such as mice and humans, where there are very few remaining commonly expressed unannotated genes.

thankyu !! what are those nUmber of secondary and non-unique alignments, ambiguous, and no feature aligned?

And what are the causes reads genomic origin - Intergenic shows higher number? sequencing error, mapping error, unannotated genes? I see in some of the post that for human genome, it's very likely cz of unannotated gens. The annotation of a genome could never be finished, as it could never be sequenced completely. Is this right?

@raju I've updated my reply.

Thanks for the answers. I have to tell you that for "alignment" with HISAT2 I didn't use any gtf file. I used only standard genome index from HISAT2 website. I used gencode new version (gencode.v27 with both protein and non-coding genes) gtf for Stringtie tool. In the above qualimap result I didn't mention the paired-end option "-pe". So, I have rerun the qualimap again with "-pe" and also "-p strand-specific-reverse" which is RF in HISAT2. The following are the results.

Reads alignment
reads aligned (left/right) = 86,704,161 / 85,485,519
read pairs aligned  = 83,397,732
total alignments = 244,302,519
secondary alignments = 72,112,839
non-unique alignments = 22,973,098
aligned to genes  = 249
ambiguous alignments = 25
no feature assigned = 9,759
not aligned = 5,887,822; Reads genomic origin

exonic =  249 (2.49%)
intronic = 693 (6.92%)
intergenic = 9,066 (90.59%)
overlapping exon = 54 (0.54%)

Transcript coverage profile
5' bias = 0.39
3' bias = 10.04
5'-3' bias = �


text:
ENA minion 2D data doesn't actually have 2D data
I am looking at some yeast ONT data from De novo yeast genome assemblies from MinION, PacBio and MiSeq platforms, Giordano et al 2017.
I noticed that the Minion fastq files have both a forward and reverse file. The methods paper says that these reads are from the SQK-MAP006 and SQK-MAP007 protocols which both appear to be 2D. Are these just the 2D forward and reverse files?
Still being unsure, I checked if this was true by counting the number of characters in each read.
>> zcat ERR1883385_1.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/1

>> zcat ERR1883385_2.fastq.gz | head -1
@ERR1883385.1 dc5b9de2-ba00-47ab-b698-8ad470704c0a_Basecall_2D_000_template EVO107591_Sanger_FAA52213_CBS432_Phenol_3a_4517_1_ch251_file44_strand/2

>> zcat ERR1883385_1.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1675

>> zcat ERR1883385_2.fastq.gz | head -2 | awk '{ print length($0); }' -
147
1693

As you can see, the read is the same but the forward and reverse strands are different lengths.
Is there some way to just download the poretools equivalent to  2D from ENA without having to download the fast5 files and extract them myself?

Probably not. ENA doesn't have any ability to specifically flag a file as 2D, and the nanopore file formats have changed a lot, so you pretty much need to rely on the submitter to separate out the reads in a meaningful way.
Also, there wasn't really any consistent convention for what ONT used in their read headers. Looking at some of my old files they had a plain sequence ID (i.e. without anything suggesting temp/comp) in different template/complement folders in the FAST5 files for the components of a 2D read and the consensus sequence, but later changed it to _2D_template and _2D_complement as in the example you have (possibly with an intermediate _1D_template and _1D_complement format). Parsing/understanding those FAST5 files can be difficult, which is why poretools was created in the first place.
It's still a bit of a mess, with the 1D$^2$ reads ending up in a Frankenstein-like state that's partially 1D and partially 2D.


text:
Filter bam using SNP list in bed format with minimum mapping quality and base quality
I have a bam file and a bed file that defines a list of SNPs. I would like to filter the bam file to contain only those reads with a minimum mapping quality that overlap at least one SNP with a minimum base quality.
Samtools seems to almost solve this problem, but not quite:
samtools depth -b snp_list -q 20 -Q 30 example.bam allows me to count the number of alignments meeting my desired criteria, but not to produce a bam with these alignments.
samtools view -L snp_bed -q 30 example.bam allows me to filter the alignments overlapping any SNP with minimum mapping quality, but not with a minimum base quality at the SNP sites. There is no -Q option for samtools view, and the -q option is for mapping quality
Can this be achieved in samtools? Or is there another preexisting tool that can?

Thanks Devon ;-)
using samjdk : http://lindenb.github.io/jvarkit/SamJdk.html
a bed file 'input.bed'
and the following snippet:
private IntervalTreeMap<Boolean> treeMap=null;
private final int BASE_QUALITY= 20;
@Override
public Object apply(final SAMRecord record) {
if(record.getReadUnmappedFlag() || record.getCigar()==null) return null;
if(treeMap==null) {
    treeMap = new  IntervalTreeMap<Boolean>();
    final String bedStr;
    try {
        bedStr = IOUtil.slurp(new java.io.File("input.bed"));
        }
    catch(java.io.IOException err)
        {
        throw new RuntimeIOException(err);
        }
    Arrays.stream(bedStr.split("[\n]")).
        filter(L->!(L.startsWith("track") || L.isEmpty() || L.startsWith("browser"))).
        map(L->L.split("\t")).
        map(A->new Interval(A[0],Integer.parseInt(A[1])+1,Integer.parseInt(A[2]))).
        forEach(I->treeMap.put(I,true));
    }

byte quals[]= record.getBaseQualities();
int refpos= record.getStart();
for(CigarElement ce:record.getCigar()) {
    CigarOperator op= ce.getOperator();
    if(op.consumesReferenceBases())
        {
        if(op.consumesReadBases()) {
            for(int i=0;i< ce.getLength();i++)
                {
                Interval t = new Interval(record.getContig(),refpos+i,refpos+i);
                if(!treeMap.containsOverlapping(t)) continue;
                int readpos = record.getReadPositionAtReferencePosition(refpos+i);
                if(readpos==0 ) continue;
                readpos--;
                if(readpos>=quals.length)  continue;
                if(quals[readpos]>BASE_QUALITY) return true;
                }
            }

        refpos+=ce.getLength();
        }

    }
return null;
}

usage:
samtools view -bu -L input.bed input.bam |\
   java -jar dist/samjdk.jar --body -f snippet.code

You're not going to get samtools to do that, though I wouldn't be surprised if something in jvarkit can do it. Baring that, here's a bit of python.
import pysam

bed = open("foo.bed")
bam = pysam.AlignmentFile("foo.bam")
of = pysam.AlignmentFile("filtered.bam", mode="wb", template=bam)
seen = set()  # memoize reads seen to not double print them
for line in bed:
    cols = line.strip().split("\t")
    for col in bam.pileup(cols[0], int(cols[1]), int(cols[2])):
        for b in col.pileups:
            if b.alignment.mapping_quality < MIN_MAPQ:
                continue
            pos = b.query_position
            if not pos:
                continue
            if b.alignment.query_qualities[pos] < MIN_PHRED:
                continue
            if b.alignment.query_name not in seen:
                seen.add(b.alignment.query_name)
                of.write(b.alignment)
of.close()
bam.close()
bed.close()

Or something along those lines, I haven't debugged that.


text:
removing indels in plink
I have data in Plink format. I removed Indels using --snps-only command in Plink, but still some snps are as below:
rs58119544;rs5746975    0   16915280    T   A

I wonder, how i can get rid of these snps? the correct should be single snps without this (;) sign. Any idea to fix it please?

Hello, I have been trying to figure out what is plink format and it seems that it does not really exist. There are a lot of plink-flavoured formats. Could you [edit] your question and be more specific? Maybe it would be also useful if you would write why are the lines with ";" a problem. If you just want to keep lines without semicolon, you can just do grep -v ";" my_plink_file > filtered_plink_file.

Use grep to get SNPs that have ";" in the name:
>$ grep ";" myMap.txt | cut -f1 -d" " > myBadSNPs.txt

Then use plink with --exclude :
plink --file myFile
--exclude myBadSNPs.txt
--out myFileFiltered


text:
Optimal design for an RNA-seq experiment with ERCC RNA Spike-In Control Mixes
I am planning an experiment on gene expression. I have 6 conditions and 3 replicates in each condition. I would like to use ERCC spike-in mixes for quality check: lower limit of detection, dynamic range and also check that the fold-changes are accurate. For example, by knowing that spike-in sequence xyz in condition 2 should be a 4-fold higher than in condition 1, does it actually get detected as a 4-fold change in my analysis? To assess the fold-changes, I need to use 2 different mixes. Only one mix should be used per sample. The manual suggests either adding mix 1 in a reference and mix 2 in all other samples or randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples? Or can I randomise among conditions (for example 1, 3, 4 get mix 1 and 2, 5, 6 get mix 2)? So far, I am intending to do the latter, so using mix 1 in conditions 1, 3, 4 and mix 2 in conditions 2, 5, 6. This would also be cheaper because I can only buy equal volumes of mix 1 and 2. Do people have any advice on this?

I have 6 conditions and 3 replicates in each condition.

If including more biological replicates has a low experimental cost, you should be using 6 replicates per condition (and fewer reads per condition, i.e. maintaining the total number of reads) in order to pick up the majority of differentially-expressed genes. For that, I'd recommend using 3 lanes of a flow cell, with 12 dual-indexed barcodes (with the same index on each end) for each lane. Doing dual indexing will protect a bit against index switching and make the results more robust.

randomly assigning mixes 1 and 2 among samples. The latter is not very clear to me. Should I use the same mix for all replicates within a condition or should I randomise really among all the samples?

Randomly assign mixes among and within the conditions. Make sure that the mixes are equally represented (as much as possible) across the conditions, but that each condition has at least one of each mix. This can be done using something like the following in R:
numReplicates <- 7;
numEachMix <- ceiling(numReplicates / 2);
assignments <- sample(rep(c("mix1","mix2"),numEachMix))[1:numReplicates]

For working out the level of fold change this is reliably replicated (or "statistically significant"), you can split the samples into mix 1 and mix 2, then apply the same algorithms for ERCC differential expression as with the non-ERCC genes.
However, as ERCC is typically added in proportion to the amount of RNA (or, if you're lucky and have cell-sorted populations, possibly the number of cells), all it's really testing is the accuracy of sample prep.

Thank you for this very clear answer. If I randomise as you say, what happens if my pseudo-random number generator assigns the same mix to all the replicates of the same condition? Shall I just switch one over to the other mix in that case? In the analysis of the fold coverages, shall I just separate my samples into 2 groups: mix 1 and mix2 and do fold-change analysis between the two groups?

If all replicates of the same condition have the same mix, then switch one of them over -- "each condition has at least one of each mix", but it'd be better to distribute the numbers of mixes equally in each group.

so, with 3 replicates, shall I do mix 1, mix 2 and mix random?

Yes, but randomly assigned. In R, you could do this by something like the following:
c("mix1","mix2","mix1","mix2")[sample(4)][1:3]

that's quite clever, I will also set a seed with set.seed(date + 1, 2, 3, etc.)


text:
How to filter info score post-imputation?
this is a very basic question, but I cannot find it explicitly stated anywhere: When exactly should I filter the info scores after imputation?
I did imputation in Impute2 and am planning on using SNPtest to analyze, but cannot tell if I should filter the scores in Impute2, SNPtest, Gtool, QCtool, or through python.
I'm hoping to get rid of the poorly imputed results (with an info score cutoff <0.5) in a population of ~4000 (half cases, half controls). After this, I'm running an association analysis -- probably through PLINK, but I'm still figuring out the best way to go from Impute2 to PLINK.
Any help would be greatly appreciated!!

In simple GWAS setups, each SNP is analyzed independently.  In those cases you can filter out SNPs with poor INFO scores at any point.
For analyses that combine information across SNPs (for example FaSTLMM), I would recommend filtering out SNPs before running the step that combines information across SNPs.
So, if you filter right after running Impute2, you should be safe.

Thanks so much, winni2k. What is the best way to filter the SNPs? Is that usually an option when inputing files to a software, or something that should be done separately (e.g., through something like python)?

Plink is wildly popular and very fast. But do use version 2 if you want to avoid surprise allele reordering.

Depending on the size of your samples (1K, 10K, 100K?) or the size of imputation region (1 gene, 1 chromosome, whole genome?) imputation files even after zipping are usually huge.
Instead of filtering SNPs and creating yet more gazillions of subset files. I prefer to have excludeSNP.txt file list of SNPs. This list then used as input for further analysis using SNPTEST flag -exclude_snps_g. This list is part of IMPUTE2 output or could be additional list of SNPs that we wish to exclude for other reasons. In short, filter at the point of analysis not the imputated files.
If info file is missing we can run SNPTEST with -summary_stats_only flag, which gives you the info score.
Here is an old still relevant post at BioStars post-imputation QC.


text:
Why is ribosomal RNA difficult to remove even with Poly(A) selection?
In this answer (actually in a comment), it is stated that:

As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite difficult to completely remove from a sample (even with polyA selection), with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to.

Which are the most relevant issues that make Poly(A) selection fail? Looking online, Poly(A) is suggested often but I didn't find an explanation for the different ways in which it can fail.

We've found ribosomal RNA to be less of a problem with sequencing that depends on polyA, which suggests the issue might be in the library preparation, rather than the selection.
Many polyA RNA library preparation methods involve amplification, rather than selection, which means that existing transcripts that are present in very high abundance (such as rRNA) will still show up in the mixture.
But polyA-tailed RNA isn't the only RNA present in a sample. For a comprehensive look at RNA, the whole-RNA sample preparation methods can't use any form of selection to identify target molecules. In this case, it's typically recommended to use a form of ribosomal RNA removal, where a probe to ribosomal RNA is attached to magnetic beads (or a column, or similar), and used to fish out out ribosomal RNA, leaving the remainder behind. Unfortunately those probes don't work perfectly, either due to the ribosomal RNA being damaged at the probe site, or the probe not matching the ribosomal RNA molecule, or the RNA just not getting into the right physical state to bind with the probe. All these things mean incomplete removal of ribosomal RNA, ending up with a situation where the highly-abundant ribosomal RNA still pokes through in the reads.


text:
Subset smaller BAM to contain several thousand rows from multiple chromosomes
There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.)
Normally I do the following, which will subset the BAM file.bam and keep the header
samtools view -H file.bam > header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam

In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2.
I could first try grepping by individual chromosome, and then combining the two SAMs
e.g. here's complete BAM with grepped chr1 and (incorrect but complete) header
samtools view -H file.bam > header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam

but then I have two problems:
(1) I may not be grepping the alignments to chromsome 2---there may exist BAM rows which contain 'chr2' but are not alignments.
(2) I think one must manually edit the header. There's probably no way around this.
Is there an easy way, Bioinformatics SO?

If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:
samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

That will select 40% (the .4 part) of the reads (123 is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the .4 part to a sufficiently small number.
In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.

You can use SAMsift:
samsift \
   -i file.bam \
   -0 'c={"chr1":5000,"chr2":5000}' \
   -f 'c[RNAME]>0' \
   -c 'c[RNAME]-=1' \
   -m nonstop-remove

Explanation:

-i file.bam – input file
-0 'c={"chr1":5000,"chr2":5000}' – initialization (create a count-down dictionary for the chromosomes of interest)
-f 'c[RNAME]>0' – filtering criterion (is the counter for the current chromosome still >0?)
-c 'c[RNAME]-=1' – code decrementing the counter of the current chromosome (5000 -> 4999 -> ... -> 1 -> 0 -> -1 -> ...)
-m nonstop-remove – remove lines causing Python errors and don't stop (in this case, an error can be caused by accessing a non-existing counter for another chromosome, e.g., for chr3)

See the SAMsift readme for more info.

If you could add some explanation about what are those options and why they are necessary. For instance I don't understand why the argument to f is c[RNAME] and not other things

I've just added the explanation.

I'd generally recommend Devon Ryan's answer. However, if you do care about having the same number from each Chr, you could use the following python/pysam code (this will output approx 5000 from every Chr):
from pysam import AlignmentFile
from random import random

nreads = 5000

infile = AlignmentFile("mybam.bam")
outfile = AlignmentFile("outbam.bam", "wb", template=infile)

for chr in infile.reference_names:
     reads_in_chr = infile.count(chr)
     frac = float(nreads)/reads_in_chr
     count = 0
     #-- Replace this loop for the exactly 5000 first reads --#
     for read in infile.fetch(chr, multiple_iterators=True):
        if random() < frac:
            count += 1
            outfile.write(read)
     print("ouputted %i reads from %s" %(count, chr))
outfile.close()

This uses all chromosomes. If you only want to use, say chr1 and chr2 replace infile.reference_names with ["chr1","chr2"].
If you want exactly 5000 reads from each chr, but don't care if they are the first ones or not, then you could replace the inner for loop with:
for read in infile.fetch(chr, multiple_iterators=True):
    count += 1
    if count >= nreads:
        break
    outfile.write(read)

if you wanted the mates of those reads as well, you could add the following after outfile.write(read):
mate_read = infile.mate(read)
outfile.write(mate_read)

Note that this is slow.

To speed things up, infile.count(chr) can be replaced with parsing idxstats() (unless one is using CRAM files).

Thanks, I'd assumed that there would be an easy way to fetch this from the index.

"To speed things up, infile.count(chr) can be replaced with parsing idxstats() " I'm not quite following. What would the above look like?

using samjdk: http://lindenb.github.io/jvarkit/SamJdk.html
$ java -jar dist/samjdk.jar --body -e \
  'Map<String,Integer> c=new HashMap<>(); public Object apply(SAMRecord r) {int n=c.getOrDefault(r.getContig(),0);if(n>=5000) return false; c.put(r.getContig(),n+1); return r;}' \
 input.bam

It'd be nice if people down-voting answers would explain why they're doing so...

@Devon (I have seen the comment until now) I downvoted because I find that a long line of code without explaining what it does internally is usually harder to understand, the code as is would improve much (IMHO) broken in several lines and with comments. It is hard for me that I don't know java to understand what is happening

If you just want to create a truncated BAM file with a header then you can significantly simplify your original code:
(samtools view -H input.bam; samtools view input.bam | head -5000) | samtools -bo output.bam

This command avoids the intermediate file and few gratuitous intermediate command invocations, at the cost of a subshell (invoked by (…)).
As above, but with formatting:
(
    samtools view -H input.bam
    samtools view input.bam | head -5000 # (*)
) \
| samtools -bo output.bam

… this can of course be extended to filter by multiple chromosomes by replacing the line marked with (*) above by one or multiple lines that subset by chromosome name (samtools view input.bam chrx, no need for grep if you have indexed the original BAM file!).

Is it samtools -bo?


text:
5' and 3' bias in Rna-seq data
I'm working with rna-seq samples. I see 5' bias and also 3' bias in the per-base sequence content plot. From this link I see that the bias at the start of the sequences appears to be the result of biased selection of fragments from the library. And this won't effect any downstream analysis. But In all of my samples I also see bias at the end of the sequences. What can be the problem? Samples are strand-specific sequenced with poly-A selection protocol.
In the Transcript coverage profile I see 5'-3' bias for all my samples is: 0.9
Showing the per-base sequence plot for one of the sample: For all the samples I see that 5' and 3' bias.
When I looked at Adapter content this is what I see:
No samples found with any adapter contamination > 0.1%

I expect there was a sequencing problem during the last base, where some of the reagents were running low on the sequencer. This won't pose any real problem, RNAseq aligners like STAR will just soft-clip the last base or two if they're mismatches.
It's common to see a bit of bias toward the 5' or 3' ends in RNAseq, mostly due to whether poly-A selection was done or there was a bit of degradation. The primary thing to ensure is that the bias is similar across samples/groups. If that's the case then you don't have to worry that your downstream analyses are going be affected by some sort of per-group bias.

I'm using hisat2 and yes I see the bias is similar across samples. Thanks for the answer.

Ah, for hisat2 you might get a bit better results if you trim off the last 1-2 bases. This is because hisat2 does global alignments (i.e., it doesn't do soft-clipping). Or at least that was the case last time I checked.

Cool. I will consider this and will trim few bases in the 3' end. Does it need any trimming in the 5' end also?

Nope, the 5’ end should be ok.

I thought tha hisat2 had a --local option, but when I run hisat2 -h | grep local, I only get --ma <int>         match bonus (0 for --end-to-end, 2 for --local). No --local proper. The help does not look coherent.

@bli The joys of bioinformatics package documentation...

I know this is an old post, but if this plot is from after trimming I would suggest a different explanation: some trimming tools remove poly-A sequences from reads. If that's the case then any read ending with A will have that removed, this leads to a 0% A base content in the final base position (and a dropping %A in the final couple of bases).
This should have no effect on alignments or analysis.


text:
Get gene exons from UCSC browser by HGNC name
I have a set of gene names in HGNC standard and I want to get all their exons (in BED format) from the UCSC browser but it seems like the browser doesn't index genes by their HGNC name. How can I do that?

Does it have to be UCSC? This is trivial with biomart.

You could grab exons via Gencode v27 (hg38) GFF annotations and convert them to BED via BEDOPS convert2bed/gff2bed:
$ wget -qO- ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gff3.gz \
    | gunzip --stdout - \
    | awk '$3 == "exon"' \
    | convert2bed -i gff - \
    > exons.bed

Then use grep to filter exons against your list of gene names:
$ grep -w -f hgnc.txt exons.bed > exons.filtered.bed

If you need a 3-column BED, you can just use standard Unix tools:
curl ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gtf.gz \
  | gzip -dc \
  | grep -wFf list.txt \
  | awk -vFS="\t" '$3=="exon"{print $1"\t"($4-1)"\t"$5"\t"$9}' > out-3.bed

If you need a 12-column BED:
# install k8 and paftools.js
curl -L http://bit.ly/k8-024 | tar -jxf -
cp k8-0.2.4/k8-`uname -s` k8
git clone https://github.com/lh3/minimap2
# actual command line:
curl ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_27/gencode.v27.annotation.gtf.gz \
  | ./k8 minimap2/misc/paftools.js gff2bed - \
  | grep -wFf list.txt > out-12.bed


text:
Strange per sequence GC content results
I would be grateful if someone could take a quick look at these FASTQC results. This is rna-seq paired-end data. From the FASTQC manual, an unusual distribution seems to be suggestive of contamination and a shift in the curve is suggestive of a systematic bias. GC content distribution both pre-alignment and post-alignment are strange.
Samples are paired end, strand specific and % of mapped reads is above 95% for all the samples. There is no adapter content also.
What could be the problem?

What are you trying to sequence? Something more specific than "poly-A amplified RNA" would be helpful to know if this is unexpected.

In the line of gringer. What specie are you sequencing ? Is there some background about GC distribution for this type of sequencing?  Did you check with other previous sequences to have such a problem? (It could be a termal microorganism and be totally normal or not...)

Its human species. except per sequence GC content and duplication levels of every sample other plots are totally fine. there is no adapter content also. samples are strand specific sequenced with poly-A selection. % of mapped reads is above 95% for all the samples. Paired end reads % average GC content is <50% or equal to 50%. % of duplicate reads is 65-70%

I see sequence bias in 5' and 3' end also in per base sequence plot. Is this GC content distribution due to bias selection?

What kit have you used to generate your libraries ?

Ribosomal rna depletion kit

Okay, so ribosomal depletion kits normally produce total RNA rather than polyA selected RNA. Normally you would do both ribo depletion and polyA selection. If you did ribo depletion rather than polyA selection, I would suggest you might be looking a contamination with ribosomal sequences. That would explain why you still see high mapping rates.

I'm not sure how normally people are combining ribo depletion and polyA, doing one obviates doing the other. However it's worth noting that ribo depletion is hit or miss regarding how well it actually works. I've seen cases with 1-2% rRNA afterward and cases with 20-30% rRNA, it depends a lot on how fresh the samples are.

@IanSudbery so is this common with ribo depletion kit? There is no adapter content also. Can I go further for with these samples?

Any kit can work inefficiently. You could check this by looking at how many reads map to ribosomal genes or ribosomal RNAs. Its not neccesasrily a problem, other than that you might end up with low counts in none ribosomal genes. I also be very careful when it came to normalisation in a DE experiment: large amounts of sequencing estate taken by a small number of genes can throw off normalisation for other genes.

Yes when I got the counts from the bam files after aligning to genome, I see most of the samples having approx. 25k-40k genes with 0 counts. Is this due to the reason what you said? How can I check reads mapping ribosomal genes or rRNAs?

@IanSudbery Hi, Could you please check my above comment. And could you please tell how to check reads mapping to ribosomal genes or rRNAs?

Sorry, our mapping pipeline does all this automatically and I had to dig out the code to remind myself how its done.
If you think you've got a very skewed counts distribution, i'd start by looking at the identities of the genes with the highest counts.
Alternatively to look directly at rRNA genes, you can download their annotations from UCSC or biomart with the correct biotypes. For the ribosomal proteins use go term GO:0003735. Then quantify with whatever you quantify with.
% reads mapping to protein coding genes is generally a pretty good QC metric.

Although generally I recommend people to be wary about the warnings that FastQC gives (it tends to be overly paranoid), the GC content graphs here do look odd. It's good that you've had a look at it with another program to confirm the observation. They should have a distribution that is fairly close to normal, a bit like this:

Without more information about your sample, I can't tell whether or not this is expected from your input. This looks like fairly substantial contamination, but could also be caused by biased selection during sample preparation (e.g. low quality sample, low yield, PCR amplification bias), or just looking at a multi-organism sample rather than a single-organism sample.
As a contamination check, I'd recommend running your RNASeq reads through a kmer-based classification tool (e.g. OneCodex, Centrifuge, BBSketch, or Sourmash). This should tell you the most likely place that contaminant sequences came from.


text:
How to create artificial or dummy phenotype data
I have whole genome sequence (human) data (SNP data), but don't have phenotype data ready for any traits yet. Hopefully, I will get real phenotype data soon. To run GWAS (genome-wide-association study) or any other analysis, i want to create artificial/dummy data for any trait, let's say human height. Could anyone please recommend what is best option? Does plink generate any dummy phenotype data? or any other programme?

You just need to create some random values (either integer for height, or categorical for "disease" or any other trait you want to "test"). What have you tried so far? Why do you need to create it now? If for learning process you might be able to download some GWAS data (take into considerations that it is identifying information so you might find restrictions to download data).

Using wakefield package it is pretty easy, here is an example from GitHub page. This will create 500 dummy variables for id, race, age, etc. :
library(wakefield)

myDummyPheno <- r_data_frame(
  n = 500,
  id,
  race,
  age,
  sex,
  hour,
  iq,
  height,
  died
)

# result
myDummyPheno

# # A tibble: 500 x 8
#    ID    Race        Age Sex    Hour           IQ Height Died
#   <chr>  <fct>     <int> <fct>  <S3: times> <dbl>  <dbl> <lgl>
#  1 001   Bi-Racial    27 Male   00:00:00     93.0   71.0 F
#  2 002   White        28 Female 00:00:00    104     64.0 T
#  3 003   Black        28 Male   00:00:00    111     67.0 T
#  4 004   Hispanic     35 Male   00:00:00     90.0   67.0 T
#  5 005   Hispanic     21 Female 00:00:00     94.0   70.0 T
#  6 006   Black        35 Male   00:00:00     74.0   73.0 F
#  7 007   Black        27 Male   00:00:00    102     68.0 T
#  8 008   White        28 Female 00:00:00     95.0   70.0 F
#  9 009   White        30 Female 00:00:00     80.0   71.0 F
# 10 010   White        30 Female 00:00:00    114     66.0 T

Awesome, that is a good pointer!


text:
seqret: Warning: bad /protein_id value
I am using the seqret tool from the Emboss suite to transform produce an EMBL annotation file from a GTF annotation file and the corresponding FASTA sequence:
seqret -sequence ecoli.fa -feature -fformat gff -fopenfile ecoli.gtf -osformat embl -auto -outseq ecoli.embl

It produces my desired EMBL file but it warns against a 'bad protein_id':
...
Warning: bad /protein_id value '"CCE57732"'
Warning: bad /protein_id value '"CCE57733"'
Warning: bad /protein_id value '"CCE57734"'
Warning: bad /protein_id value '"CCE57735"'
Warning: bad /protein_id value '"CCE57754"'

What's wrong with these protein IDs? The organism I am working with is Escherichia coli.
Here are a few lines from my GTF file. The first line has a protein_id.
Chromosome      ena     CDS     349     2808    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding"; protein_id "CCQ26880";
Chromosome      ena     start_codon     349     351     .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";
Chromosome      ena     stop_codon      2809    2811    .       +       0       gene_id "HUS2011_0001"; transcript_id "CCQ26880"; exon_number "1"; gene_name "thrA"; gene_source "ena"; gene_biotype "protein_coding"; transcript_name "thrA-1"; transcript_source "ena"; transcript_biotype "protein_coding";

To me it looks like a bug, where the encoding or the quotes are not properly handled. (ie, CCE57754 would be a valid protien_id(?) but "CCE57754" it isn't)

Can you show (a small part of) your GTF file? I tend to agree with Llopis, but not all fields in the GTF files can be quoted.

You're supplying a GTF feature file, but telling it that the features are in GFF3 format. I think you need -fformat gff2. See https://www.ebi.ac.uk/Tools/sfc/emboss_seqret/help/

I tried with gff2 but I get the same problem: Warning: bad /protein_id value 'CCQ26880'
Warning: bad /protein_id value 'CCQ26881'
Warning: bad /protein_id value 'CCQ26882'
Warning: bad /protein_id value 'CCQ26883'

If your copying the output directly, that's a different problem ('CCQ26880' vs. '''CCQ26880"'). It looks like it's parsing the GTF correctly but choking on the values for some reason. I'm not very familiar with this tool so I can't  advise further.

Do you use another tool to get EMBL files from GTF and FASTA?

@charlesdarwin have you tried this to convert FASTA -> EMBL? http://sequenceconversion.bugaco.com/converter/biology/sequences/embl_to_fasta.php (no idea of how it works, but you can try)