text:
How to create a .bed file from .fasta?
I have some problems in creating a .bed file for hg19, so I will be able to visualize the .bed file in IGV.
The .fasta file contains rows of this form:
>chr1:0-1000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:1000-2000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:2000-3000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>chr1:3000-4000
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

I tried both with biopython:
from Bio import SeqIO

with open(f1) as in_f, open(f1+'.seq','w') as out_f:
    for record in SeqIO.parse(in_f, 'fasta'):
        out_f.write('{}\t0\t{}\n'.format(record.id, len(record))

and awk:
samtools faidx file.fasta
awk 'BEGIN {FS="\t"}; {print $1 FS "0" FS $2}' file.fasta.fai  > file.bed

Both ways gave me the same output which IGV couldn't work with:
chr1:0-1000 0   1000
chr1:1000-2000  0   1000
chr1:2000-3000  0   1000
chr1:3000-4000  0   1000
chr1:4000-5000  0   1000
chr1:5000-6000  0   1000
chr1:6000-7000  0   1000
chr1:7000-8000  0   1000

Now I tried to see what would be a valid .bed file format from UCSC and the format looks different. For example:
#bin    chrom   chromStart  chromEnd
13  chr1    39845715    39846245
19  chr1    87030756    87031903
22  chr1    111148649   111149175
22  chr1    113246065   113246628

Why my methods don't output the .bed in the same format?
How can I make the awk/biopython methods to output a .bed file that IGV can actually load.

As @devon-ryan points out below, both of your BED format examples are incorrect. I think it would be helpful for others if you edited your question to give a correct example. See for example http://genome.ucsc.edu/FAQ/FAQformat#format1

Here is a fun little python script I cooked up for the occasion. It will take a standard fasta file as a command-line argument and turn it into the proper bed format. Using your example fasta:
import pandas as pd
import sys

inFasta = sys.argv[1] # take fasta as command argument

def fastaParser(fasta):
    headers = []
    with open(fasta) as f:
        header = None
        for line in f:
            if line.startswith('>'): # identifies fasta header line
                headers.append(line[1:-1]) # append all of the line that isnt >
                header = line[1:] # in reset header
    newHeader = (header.replace(':',',') for header in headers) # format to be accepted later
    newnewHeader = (header.replace('-',',') for header in newHeader) # format to accept later
    bedHead = (header.split(',') for header in newnewHeader) # separate by comma from format above
    return bedHead

fastaParsed = fastaParser(inFasta) # run the function

headers = list(fastaParsed) # go from generator to list
bedfile = pd.DataFrame(headers) # create dataframe which will output bedfile

bedfile.to_csv('fastabedfile.bed', sep='\t', index=False, header=None) # output bedfile format

to run:
[dkennetz@nodecn202  fasta_parser]$python fastaParser.py fakeFasta.fa

[dkennetz@nodecn202  fasta_parser]$more fastabedfile.bed
chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

I'm assuming the FASTA header contains all the information you want. You can easily generate a proper bed file with Biopython:
from Bio import SeqIO

for record in SeqIO.parse('test.fa', 'fasta'):
    lst = record.id.split(':')
    chr_, (start, end) = lst[0], lst[1].split('-')
    print('\t'.join((chr_, start, end)))

Note that native Python code is only marginally longer:
with open('test.fa') as f:
    for line in f:
        if line.startswith('>'):
            lst = line.rstrip().split(':')
            chr_ = lst[0].lstrip('>')
            start, end = lst[1].split('-')
            print('\t'.join((chr_, start, end)))

Both code snippets output bed format:
chr1    0   1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

The UCSC format you linked to isn't a BED file, your method should never produce it. What you posted as your desired output is also not a BED file. Below is a BED file:
chr1    0       1000
chr1    1000    2000
chr1    2000    3000
chr1    3000    4000

awk and biopython were producing the appropriate thing already. You just need to zoom in on a region in IGV. Since it's not a bigBed file, IGV won't show everything at all zoom levels.
As a demonstration, below is the BED file from above in IGV:

If your headers are always of the form >chrZ:a-b, where chrZ is a chromosome name, a is a zero-indexed, half-open start position, and b is a zero-indexed, half-open stop position:
$ awk '($0 ~ /^>/){ print substr($0,2) }' in.fa | sed 's/[:|-]/\t/g' > out.bed

If your headers contain coordinates using a different index, then you can make adjustments to $1 and $2 in the awk statement to create zero-indexed, half-open coordinates to export to BED, depending on the initial index.

How do I do it with biopython?

I would not use Biopython for this. I recommend using standard Unix tools that are fast and flexible. I would not use a hammer to open a jar.

Assuming that you have a reference genome sequence file (e.g., reference.fasta) available. Getting this kind of file is straightforward.
1) Index the reference genome and map your reads or FASTA sequences to it (for example with bowtie2)
# index reference genome (should be precomputed)
bowtie2-build reference.fasta reference
# map reads
bowtie2 -x  reference -U file.fasta -S file.sam
# compress SAM to a BAM (binary) file
samtools view -bS file.sam > file.bam

2) Create a BED files from the BAM file using bedtools
# extract BED file
bedtools bamtobed -i file.bam > file.bed

Alternatively, for sliding windows you can generate these from a reference sequence provided that you know the length of each chromosome (perhaps there is a way to extract these directly from reference.fasta):
# length per chromosome
samtools view -H file.bam | grep "SQ" | cut -d":" -f2-3 | sed 's/LN://' > file.chr.txt

Create sliding windows (using only the lengths per chromosome)
# create sliding windows
bedtools makewindows -g file.chr.txt -w 1000 > windows.bed


text:
Compare two networks of the same genes between two species
I had a set of genes for honey bees. I found the homologues of those genes in the house fly and created a network of the honey bee genes and a network of the homologous fly genes using STRING.
I know that I could do a topological comparison between the networks and check if the same set of genes give the same hubs when in two different species. However, I can't find a single study where something similar has been done.
What biological inference can I draw from the topological parameters? For example, if the diameter of the fly network is smaller than the honey bee network, what inference could I draw?
UPDATE: It would be great if atleast some publications could be given where network comparisons have been done between different species.

So, you can give a try to this Cytoscape app that will allow to perform comparison between networks.
This app is not specifically designed for comparing networks across different species but it will enable you comparing any kind of networks.
These are the papers they published, which can give you some hints.
About your specific problem...I can't help you about what to look for when you're comparing networks from different species (not an expert in this specific topic) but as far I know you should look for something related to patterns occurring across different networks (i.e. species), network motifs, distribution of network centralities, conserved structures. In general, to compare networks, I would look for global properties (e.g. diameter, clustering coefficient...) while I would use local properties (e.g. degree, closeness, and network centralities in general) to compare the role of a node that you retrieved in all your networks.
Also, since you're looking for some literature, these links may help:

About conserved patterns of protein interaction in multiple species
About comparing genomic expression patterns across species
This one is slightly out of topic but can give you
some other ideas about comparing networks.

source: Stack Exchange

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

sort: votes


text:
Error compiling HMMER v3.2.1
I tried to compile HMMER v3.2.1 in Cygwin using
./configure

followed by
make

or
make install

but failed. I got this config.log file to troubleshoot and I suspect the solution is in this line:
configure:2908: error: no acceptable C compiler found in $PATH

What does it mean? Is it really all I need to address? I followed the installation steps in this document.

It means you need to download a C compiler. I would recommend using gcc!

Have a look here on how to install gcc https://askubuntu.com/questions/271388/how-to-install-gcc-4-8

gcc is a C compiler

I would strongly recommend to use Bioconda instead of compiling. It is one of the channels ("distributors of packages") in Anaconda Cloud which is basically a package source for conda.
In short - conda allows you to create virtual environments isolated from system and to install various precompiled software packages with one command (so, by the way, you can easily use multiple software versions on one system). Also, Bioconda is just one channel aiming on various bioinformatics software, but many more packages can be found in aforementioned Anaconda Cloud.
Tutorial how to start using Bioconda is on its homepage, and it is really easy (far more easier than to struggle with compiling).
Bioconda HMMER package
P.S. Of course there are cases in which you need a specific build of given software. Than it is really needed to compile manually.

Thank you. The website says "Bioconda supports only 64-bit Linux and Mac OSX." I am a Windows user, so I assume I will need to go back to compiling. Bioconda sounds amazing though; thanks for bringing it to my/our attention.

I tried to set up channels with conda config --add channels defaults, conda config --add channels bioconda & conda config --add channels conda-forge in my Anaconda prompt which seemingly worked, but conda install hmmer didn't.

@JernejT On Windows? And what does conda install hmmer output? Anyway, Linux is much more supported platform for bioinformatics. Worth to try :)

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


text:
How to paste RSIDs in CADD output
I want to paste RSIDs in CADD output as CADD does not give RSID column in its output. For this purpose I am using bedtools intersect to compare two files and have RSID column in my CADD file.
This is the command I am using.
bedtools intersect -a CADD.bed -b Input.vcf| cut -f 3

But this command is giving me error
ERROR: file CADD.bed has non positional records, which are only valid for the groupBy tool.

Is there any way to solve this error or any other way to get RSID column in CADD output file.

You might consider vcfanno.

It would be useful if you show us some line of your CADD.bed and Input.vcf and how your desired output should look like.
bedtools is complaining about non positional records in CADD.bed, which means there are lines in this file where the first three columns aren't CHROM, START, END.
I guess you will have some header lines in the beginning that aren't introduced by an #. Just delete them (or add in # to the beginning of the line) and try again.
fin swimmer


text:
Is there a read mapper that allows input in SAM format instead of FASTA/FASTQ?
When we get raw sequenced reads in FASTQ format, we usually need to do some pre-processing (QC, demultiplexing, etc.) prior to sequence alignment. The only way to register the results of this preprocessing is to add metadata to each read. This is best done by using the SAM format which has a well-defined system of metadata tags.
However, it seems that not a single read mapper (except maybe for bwa aln with option -b) accepts SAM (or one of the binary equivalents BAM/CRAM) as input to make use of the metadata that has been added during preprocessing. This implies that one has to convert the unaligned reads back to FASTQ format.
In any case, whether one uses FASTQ throughout or converts back to FASTQ, a common annoyance with FASTQ is that metadata is registered as part of the query name - this is not standardized or even well-defined. After alignment, we again have to do some text processing to get the metadata out of the query name and into proper metadata tags defined by SAM. This back and forth hacking is highly unfortunate and error prone. I have seen this, for example, in some single-cell sequencing projects where we have highly multiplexed sequencing runs with reads that need to be QC'd and classified. The new SAM standard from May 2018 does indeed define proper tags for handling cell-multiplexed runs. SAM is the de-facto standard to store both unaligned and aligned reads for years to come so it should be standard for read mappers to accept it as input.
I would be grateful for some hints and if we could collect here some read mappers that do accept SAM/BAM/CRAM as input format.

Why not store the read identifiers and metadata in a separate table file?

Of course, this is possible. But apart from not being an answer to my question, this is not a good solution because it is not standardized and invites inconsistencies. You have to write special code to read the table. The SAM format (which is also a table) is maybe not the ideal but it is the best we have. If we could all stick to the SAM format instead of introducing custom formats the interoperability of bioinformatics tools would be much better.

I’m a bit confused because the preprocessing you mention is pretty routine, and is generally done without the involvement of BAM files.

Please remember to confirm an answer once you've received one.

The common solution for scRNA-seq is to put cell barcodes and such in read headers and then post-process things with UMItools.
But regarding your actual question, STAR can accept SAM/BAM as input with the --readFilesType SAM PE option (it's SAM for both that format and BAM). You can swap SE for PE if your data is single-end (or is effectively that way, as is the case in many scRNA-seq experiments). Note that in the case of PE experiments that mates must be next to each other. I really don't know what this does with auxiliary tags though, so you would need to see if they're lost before you start processing a lot of data).

With bwa-mem or minimap2, the recommended way is
samtools fastq -T BC,RX name-grouped.sam | bwa mem -C -p -

This passes the BC and RX tags through bwa-mem and copies them to the SAM output. If you are running paired-end alignment and your BAM is coordinate sorted, you have to run samtools collate to group reads by name:
samtools collate -Ol0 in.bam tmp | samtools fastq -T BC,RX - | bwa mem -pC -

Here bwa mem -C or minimap2 -y copies FASTQ comments to SAM output. This strategy is much easier to deal with than SAM/BAM. A caveat is that it doesn't work well for read groups because each read group requires a separate header line. You have to manually modify the header.
I am not sure which short-read RNA-seq mappers have implemented a similar feature. You might consider to send a feature request.

At the Broad Institute, we do something similar to what is described by user172818. We use uBAM (unmapped BAM) to store unmapped reads, because as you say the SAM format has a more appropriate way to encode metadata. However, we use a Picard tool called SamToFastq to do the initial conversion then pipe the output of that to BWA mem. Then we use an additional step to reattach metadata, read group info and so on, using a Picard tool called MergeBamAlignments.
You can see the full implementation here (starts as line 71):
https://github.com/gatk-workflows/gatk4-genome-processing-pipeline/blob/master/tasks/Alignment.wdl

SNAP is a hash-index based genomic aligner that can take SAM and BAM as input. It has been used as an aligner for RNA-seq (see SNAPR) and BS-seq (see BS-Seeker3). SNAPR seems to accept SAM/BAM but it is not clear whether BS-Seeker3 wraps this functionality as well.
I tried SNAP and it seems to do the job well, i.e. it only adds tags during alignment and does not remove metadata that has been deposited prior to alignment, neither from the SAM header nor from individual reads.
NOTE I post this answer to my own question after I became aware of SNAP.


text:
Downloading multiple SRA files from several SRA accession IDs does not work
I am trying to download multiple SRA files located in several SRA accessions. Some of my accession numbers are as follows:
SRR1002678
SRR1158898
SRR1158907
SRR1158950
SRR1181152
SRR1181155
SRR1181196
SRR1181204
SRR1181205
SRR1181208
SRR1269497
SRR2100227
SRR2100419
SRR2100830
SRR2101344
SRR671725
SRR671726
SRR671740
SRR671759
SRR671760
SRR671761
SRR671779
SRR671782
SRR671815
SRR671821
SRR671842
SRR671862
SRR671863
SRR671864
SRR671865
.............

I tried the following command:
for i in {SRR1002678..SRR1184123}; do wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR1*/$i; done

but this did not work.
I tried prefetch as well:
sort -u SRR.txt | parallel -j 4 "prefetch {}"

where SRR.txt is the list of my SRA IDs. However, with this command only a few files are downloaded as *.sra files, while with several other SRA IDs I get the following error:
2018-11-15T06:17:20 prefetch.2.8.2: 1) Downloading 'SRR1183962'...
2018-11-15T06:17:20 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:21 prefetch.2.8.2: 1) 'SRR1183962' was downloaded successfully
2018-11-15T06:20:27 prefetch.2.8.2: 'SRR1183962' has 0 unresolved dependencies
2018-11-15T06:20:38 prefetch.2.8.2: 'SRR1183962' has no remote vdbcache

2018-11-15T06:18:59 prefetch.2.8.2: 1) Downloading 'SRR1184034'...
2018-11-15T06:18:59 prefetch.2.8.2:  Downloading via https...
2018-11-15T06:20:40 prefetch.2.8.2: 1) failed to download SRR1184034
2018-11-15T06:20:40 prefetch.2.8.2 sys: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module - mbedtls_ssl_read returned -76 ( NET - Reading information from the socket failed )
2018-11-15T06:20:40 prefetch.2.8.2 int: libs/kns/unix/syssock.c:606:KSocketTimedRead: timeout exhausted while reading file within network system module -

I do not have a range of SRA files to download but quite random SRA files to download. What am I doing wrong? Any help is highly appreciated.

First of all, the brace expansion you tried to use isn't valid. While bash can expand {1..3} to 1 2 3, it has no way of knowing how to expand {SRR1002678..SRR1184123}:
$ echo {1..3}
1 2 3
$ echo {SRR1002678..SRR1184123}
{SRR1002678..SRR1184123}

Perhaps you meant SRR{1002678..1184123}? But there's no need for this, you already have the accessions in a file and can read that. I see that the FTP server you are accessing has organized the data like this:
sra/sra-instant/reads/ByRun/sra/SRR/SRRNNN/SRRNNNXXXX/SRRNNNXXXX.sra

Where NNN are the 1st three digits of the accession. So, for example, the file for accession SRR1002678 is at:
ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/SRR118/SRR1181152/SRR1181152.sra

With this in mind, you can write a little script that extracts the URL from each accession:
while read i; do
    n=$(echo "$i" | sed -E 's/(SRR...).*/\1/');
    wget ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/SRR/$n/$i/$i.sra
done < SRR.txt

The only trick is the sed command which simply looks for SRR and the three next characters and removes everything else on the line (it "captures" the found pattern and then replaces the line with what was captured (\1)).

Thanks @terdon but im getting the same kind of error please see:
SRR1158898%0D/SRR1158898%0D.sra            => ‘SRR1158898%0D.sra’ Resolving ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)... 130.14.250.11, 2607:f220:41e:250::13 Connecting to ftp.ncbi.nlm.nih.gov (ftp.ncbi.nlm.nih.gov)|130.14.250.11|:21... connected. Logging in as anonymous ... Logged in! ==> SYST ... done.    ==> PWD ... done. ==> TYPE I ... done.  ==> CWD (1) /sra/sra-instant/reads/ByRun/sra/SRR/SRR115/SRR1158898\r ... No such directory ‘sra/sra-instant/reads/ByRun/sra/SRR/SRR115/SRR1158898\r’.

@SafinaA.R looks like your file was created in Windows, not Linux. Note the \r at the end of the path in the error message. Fix the line endings by running sed -i 's/\r//' SRR.txt and try again.

@SafinaA.R please don't do this in the comments. They are very hard to read and shouldn't be used for discussion. Either edit your question and show the commands you tried and how they failed or come into our chat room.

For me it works best to combine cat, while and official download tool fasterq-dump.
cat your_accession_list.txt | while read i
do
fasterq-dump ${i} <options>
done

As commented above, your_accession_list.txt should have Linux endings. Link to fasterq-dump is here.


text:
Calculating alignment/mapping time
I  am trying to assemble a plant genome using AWS resources using velvet. Plant genome is huge (> 10 times human genome) and coverage is around 30 x. We are planning for de novo assembly with Velvet (threads enabled).
I would like to know if there is a calculator that can provide approximate time taken for assembly, given the necessary details. For eg. if I furnish RAM, CPU (type of instance), genome size, approximate coverage, type of sequencing (PE or SE) and number of client nodes, it should give me the number of hours or days that would take for assembly (reference and/or de novo).
Note: Previously posted on biostars, but no replies.

I don't know much about assemblies, but do you know how complex it is? chromosomes, repeated regions and other information about the complexity of the DNA? Are you using long reads or short reads (I think that short-reads, but just to make sure). Lastly, could you provide the link to the post on biostars? (Just so if someone answers there people who find this post know where to look for)

BioStars link: https://www.biostars.org/p/348943/

IMO this is unknowable, which is why you received no reply on biostars.

Velvet won't work for a 30Gb genome for sure. You might have a better chance with assemblers capable of human assembly, such as soapdenovo and discovar. But even with them, probably you won't get reasonable results.


text:
Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
After alignment using BWA, I have removed the dupliment using the samtools(Version: 1.9).
My procedure is as follows:
bwa mem -k 32 -M ref.fa read1 read2 > out.sam
samtools view -@ 0 -b -T ref.fa -o out.bam in.sam
samtools sort -n -o out.nameSrt.bam in.bam
samtools fixmate -r -m in.nameSrt.bam out.fixmate.bam
samtools sort -o out.fixmate.sort.bam in.fixmate.bam
samtools markdup -r -S -s in.fixmate.sort.bam out.markdup.bam
samtools flagstat in.markdup.bam > out.markdup.flagstat

The flagstat output result is as follows:
21611397 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
21611397 + 0 mapped (100.00% : N/A)
21422330 + 0 paired in sequencing
10711165 + 0 read1
10711165 + 0 read2
19797684 + 0 properly paired (92.42% : N/A)
21422330 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
1306000 + 0 with mate mapped to a different chr
727043 + 0 with mate mapped to a different chr (mapQ>=5)

Why is the total read number still more than the paired in sequencing after removing the duplicate in samtools flagstat output?
Is there anything wrong with my procedure?

Bwa-mem may produce chimeric alignment, where different parts of a read are mapped to distinct loci. flagstat counts them as two reads.


text:
Exporting nt's differences from MEGA alignment
I have an alignment of some full genome sequences of DENV1 that we have recently sequenced from the serum of patients. I have already cleaned the sequences, built the alignments and done the SNPs calling. After that I got the consensus sequences from each sample and I used MEGA to build the alignments of the consensus sequences and to construct a tree etc.
From just looking the alignment of the consensus sequences I can see several nt s differences among the sequences.
Does anyone know if you can export from MEGA like a sort of table where you could see the nt changes for each position and sample?

You don't say which version of MEGA you are using, but this works for MEGA X on Linux:

Open your alignment with the "Sequence data explorer"
At the "Highlight" menu, select "Variable sites".
Then at the menu "Statistics", select "Use only highlighted sites" and select what you consider an appropriate output - I choose "Display results as CSV".
Now at the menu "Data" select "Export data"; at the "Sequence data export options window" if necessary select the data format again, and select "For each site" at "Writing site numbers", and select "Only highlighted sites" at "Selected sites to include".

This will result in a comma-separated values file, showing only positions with variable sites.


text:
How to select point mutations from MAF file
I am using a mutation MAF file from GDC.
I want to select for point mutations, therefore I selected Variant_Classification equal to SNP, and Variant_Classification equal to INS and DEL - but in this case I need to check against the size of the indel mutation.
How can I check for point mutations in the indel case? I think I have to work it out through checking the length of the fields in columns Reference_Allele, Tumor_Seq_Allele1, Tumor_Seq_Allele2 but I am not sure if this is the best solution.


text:
Finding orthologues using BLAST on the NCBI database
I'm an informatics student who has essentially zero knowledge of biology. I BLASTed my gene and have 1000s of results with very low E values. Where do I go from here if I want to find orthologues?

Please [edit] your question and give us some more details. What type of blast? Blastn? Tblastn? Blastx? Something else? Is your database protein or nucleotide? And your query? What species are you working with? What kind of orthologs do you want? (see here for a short discussion of orthologs vs paralogs; I suspect you just want general homologs).

BLAST is often a sensible way to find gene/protein homologs, but resolving those as orthologs vs paralogs etc. is a non-trivial task, especially when considering large numbers of species.
Since you are new to bioinformatics I recommend you search an ortholog database that already exists.
Here are some databases:

orthoDB
OMA
EggNog

All of these should allow you to search either with text or a sequence to identify the orthologous groups that your input belongs to at different taxonomic levels. Sequence searches either follow pairwise alignment (e.g. BLAST, DIAMOND, RAPSEARCH) or profile HMM type (e.g. HMMER3) approaches.
A more complete list of ortholog databases is available from "Quest for orthologs".

Two solutions based around the same premise:

Repeat the blast using "blastp". I suspect you have used blastn (nucleotides). To translate the gene Expasy is okay and find the common "translate", here.
Blasted via Blastx without any modification to the gene.

Blast is capable of doing just about anything if you are aware of all its parameters.
The reasoning is that the gene you have will be under-represented in the database. By switching to amino acids there is a much greater chance of finding orthologous hits.

Without jumping deep into the rabbit hole of orthologology, you can just use biomart
An instruction is here


text:
ANY (technical) reason behind submitting sequences to GenBank versus ENA Sequence
The DNA sequence sections of the three INSDC databases (i.e., DDBJ, ENA Sequence and GenBank) are synchronized periodically and strive to keep their stored data as ubiquitously accessible as possible. Except for idiosyncrasies in their data submission routes, there should be little, if any, reason for preferentially submitting sequence data to one database over another. Yet many researchers display very noticeable preferences!
When asked point blank, many colleagues emphasize that these preferences are not linked to the database-specific data submission or retrieval interfaces available. Yet they cannot provide a technical reason for their choice either.
Maybe researchers base their preference on geographic vicinity ("buy local"), but maybe there are genuine technical reasons for their preference (even if they are unaware of these reasons themselves).
Thus, my question is:
Can you think of any technical reason - however minute or seemingly insignificant - in preferentially submitting DNA sequences to NCBI's GenBank over ENA Sequence (or vice versa)? For example, there may be differences in data storage or accessibility that are relevant to you.
Put differently: Short of flipping a coin, why would you (as a bioinformatics-prone end user) select one database over the other for your data submission?
--
Edit 1: Many SE posters attempt to reassure me that I can submit my DNA sequences to either database and that "the data will be just fine there". Having submitted data to both GenBank and ENA Sequence for 10+ years myself, this is not what this question is about. Instead, this question is about carving out genuine technical differences from a user perspective.
Edit 2: Since no answer has been posted yet despite a bounty, I figure this is a hard question as specified. Thus, I will broaden the question and now also accept non-technical reasons for preferentially submitting DNA sequence data to either GenBank or ENA Sequence. (That being said, I still hope for someone listing a technical reason.)

@David I am deliberately posting this question at SE Biology because I wish to query the biologists' opinion on this. The bioinformatics opinion I can provide myself. Can you think of any difference in data availability or retrieval, however small or insignificant?

The problem is that questions the answers to which are a matter of opinion are off topic here. But, a comment. I last submitted sequences over 35 years ago and at the time I used EMBL (now presumably ENA) as it was European. However as nowadays I and most others only use NCBI to look for sequences, I really see no point submitting anywhere else. I have recently had to submit RNAseq data to ENA (under a grant obligation) and the user interface for that I find truly dreadful. The Yanks have more money, so all along NCBI has had better user interfaces than EBI. Take advantage of it.

@David Thanks for your response. Actually, your reference to a recent submission to ENA due to "a grant obligation" is quite relevant to me. Having certain restrictions on the data submission imposed by a funding agency can be a very important, albeit non-biological, reason for submitting specifically to ENA. This is the type of info I am looking for. Can you provide a bit more details on this (e.g., is this restriction standard practice for some funding agencies)?

Probably best to take this off list. If you check my bio and personal website you can easily work out my University affiliation and get my email from their website.

By your comments you seem to ask about differences between managing the data repositories. Perhaps you could [edit] the question to change the title and refocus the question

@llrs I refocused the question.

@David and Michael, you can also take this to chat here, if you find that easier. Either in our site's main room or by creating a separate room for it. Just go to this site and click "create new room".

note that NCBI has both genbank and SRA (sequence read archive). SRA is for next-gen sequencing, and is a different database and submission than genbank. EMBL similarly has ERA. They're also cross referenced. Word on the street is that ERA is more user-friendly for data upload, and has less strict metadata requirements.

Opinion: I prefer to use EBI/ENA because the turnaround time for "customer service" on editing metadata about a submission is faster than on NCBI. Inquiries to change metadata on NCBI can go unanswered for up to six months in my experience.
Answer to your question: No, there are no technical reasons that make it preferable to submit data through one website or the other. All entries end up being hosted by both databases anyway. The preference is up to the user to select whichever upload interface they prefer.


text:
calling diploid SNVs from long reads
I'd like to call diploid SNV variants from long-reads data (~80SMRTcells PacBio).
I have generated a draft reference genome for an indivudual from a heterozygous (~4%) species (Canu+Haplomerger2).
I can use this reference for some applications:

RNAseq mapping;
functional annotation;
duplicate detection;
synteny analysis

but not others:

analysis of compound heterozygotes;
allele-specific expression;
variant linkage

This is mainly why I need phasing. FYI, I tried Falcon + FalconUnzip too, but the daligner phase was just taking too much time/resources and I need to find a strategy to make it work in reasonable time (working on it).
Meanwhile, for phasing, It was suggested to use HapCUT2 - from Canu's FAQ

(1) "Avoid collapsing the genome" in Canu to get raw diploid assembly as complete & separate as possible (2) HaploMerger2 to get high-quality reference and alternative haploid assemblies (3) HapCUT2 or other phasing tools to get the high-quality haplotype assembly based on the reference haploid assembly.

However, HapCUT2 requires the input of a diploid SNV file. Which could be done if I had short-reads data, but yet haven't found a method for doing so with in my hands just long-reads data.
Strategies I tried to get the SNVs:

marginPhase : method paper

but I had a issue in running it

Clairvoyante : method paper

but I currently don't have access to suitable GPU machines to build the model - and I can't use the models provided by the authors because my organism is too different from human

Any input would be greatly appreciated!

At 4% heterozygosity, falcon-unzip is probably the best hope. Trio binning would be better but you don't have the right data.

@user172818 is it possible in your knowledge to run Falcon-Unzip on non-Falcon output? From here it looks like it's a no.

You may want to try the longshot tool (https://github.com/pjedge/longshot) developed for calling variants in diploid genomes from long read data.

You could try aligning your reads to the draft reference genome with for example minimap2 and calling variants with freebayes. It appears there is a protocol for long reads: https://github.com/ekg/freebayes#re-genotyping-known-variants-and-calling-long-haplotypes
Disclaimer: I have not run freebayes myself, and it may be that this is a terrible idea. I would expect that calling SNVs (Single-nucleotide variants) from PacBio reads is very challenging.  However, with enough coverage, I think it could be possible.


text:
full visualisation of draft genomes alignment
I have two draft genomes (aprox. 500Mb each), one with ~10'000 scaffolds (genome1.fasta) and one with ~1'000 contigs (genome2.fasta). I would like to visualize their pairwise alignment, which a while ago I found it could easily be done by mummerplot in MUMmer.
In mummer, pairwise draft genomes alignment is supported:

Many times it is necessary to align two genomes that have not yet been completed, or two genomes with multiple chromosomes. This can make things a little more complicated, since a separate alignment would have to be generated for each possible pairing of the sequences. However, both NUCmer and PROmer automate this process and accept multi-FastA inputs, thus simplifying the process of aligning two sets of contigs, scaffolds or chromosomes.

It is still ok if I do it aligning small portions (e.g. contig vs scaffold), as I could get here. However, if I do it all-vs-all, I end up with a "confetti"-plot (self-coined term) instead of a dot-plot :

I am not sure how to interpret this, and what I would like to see is more of a straight line as the one. Probably it's just a matter of filtering, so I'd like to ask for advice.
I used nucmer through delta-filter as follows:
# mummer version //4.0.0beta1

# nucmer
nucmer --prefix=ref_qry genome1.fasta genome2.fasta

# ref_qry.delta generated from previous step
delta-filter -r -q ref_qry.delta > ref_qry.filter

# plot
mummerplot ref_qry.filter -R genome1.fasta -Q genome2.fasta --filter --layout

Reasons for options -r -q in delta-filter step :

To determine a one-to-one mapping of each reference and query, combine the options and use delta-filter -r -q.

I'll also accept other methods for visualisation (doesn't need to be withh MUMmer). Any input is greatly appreciated!

EDIT
I found some solutions which I reported as answers. I will accept anyone giving answer on how to tune MUMmer parameters to achieve the goal I describe in the original question.

http://dgenies.toulouse.inra.fr/

Found a solution, using D-Genies, worked great.
Some examples from their website:

Thanks to @user172818.

Also, Mashmap has an embedded visualisation tool generating dot-plots to visualize mappings, ideal for:

highly similar species or to compare different assemblies of the same species.
mapping genome assembly or long reads (PacBio/ONT) to reference genome

An example from their website:


text:
What are fundamentals to work on bioinformatic data analysis ?
Well, I'm a MSc student in Information Technology Engineering, trying to become more familiar in the Complex networks specially Bioinformatic field to get more knowledge of cancer detection by Biochemical data analysis.
beside python / R language what else should I work on to prepare myself for an academic seminar?

You need to be much more specific in order to get helpful answers. What analysis are you planning to do and which tools are you currently using? If you want a broader overview of the field, there are many courses, books, and blogs on this topic.

Which academic seminar (About which topic? what are the ponents?) ? What do you know about biology and health? Are you familiar with experiments performed on the laboratory? BTW this question is quite broad, we can answer about a concrete answerable questions. You could try asking at biostars.org but please add more information to the question in any case


text:
Question about the dots on Quartile groups in boxplot
I have Microarray Normalized Expression data for a specific Gene. It looks like below in a dataframe B
SampleID    Gene    Type
Sample1     5.02    Tumor
Sample2     5.06    Tumor
Sample3     5.1     Tumor
Sample4     5.11    Tumor
Sample5    5.127    Normal
Sample6     5.12    Normal
Sample7    5.138    Normal
Sample8    5.149    Normal

I see that the minimum expression value in the table is 5.0 and maximum expression value is around 5.9.
I wanted to show the expression between two conditions with a boxplot and used following code.
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               add = "jitter", ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

This gave a plot like this
But when I remove the jitter from the code,
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))
q + stat_compare_means(method = "t.test") +
  geom_point() +
  stat_n_text()

I see many black dots in the middle of the boxplot like this
May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?

I don't know from where does the ggboxplot function comes from or stat_compare_means? Also please, continue [edit]ing the question to clarify the question. It is still unclear to me if you want to remove some points, add them or something else

from the package ggpubr

I am confused by the wording in "May I know why I see those dots in the second boxplot after removing the jitter. Is there a way to avoid that?" Are you trying to ask  "Why did the points not disappear after I remove jitter?"

You add the points with geom_point(). Just remove it and you will get your "empty" boxplot.
q <- ggboxplot(B, x = "Type", y = "Gene",
               color = "black", palette = "npg",
               ylab = 'Gene expression', xlab=FALSE,
               order=c("Normal", "Tumor"))

Unfortunately I couldn't use stat_compare_means(method = "t.test") and     stat_n_text(), but this two just add the labels.

yes, this is the mistake I'm doing. I removed geom_point(). It looks fine now. thanq


text:
What is the "primer masterfile" format?
I am working with an amplicon kit from Swift Bioscience and they provide the primer locations in what they call a "masterfile". The format is like this:
$ head -n2 Acel-Amplicon_56g_masterfile.txt
1       43814951        43815064        MPL_S505N-A519T_F/2     43814932        43814951        MPL_S505N-A519T_F/2_primer      43815064  43815081 MPL_S505N-A519T_R/1_primer      GCCTGGATCTCCTTGGTGA     CTTCGGCTCCACCTGGT
1       115256480       115256576       NRAS_G60-S65_new_R/1    115256457       115256480       NRAS_G60-S65_new_R/1_primer     115256576 115256602        NRAS_G60-S65_new_F/2_primer     TGATGGCAAATACACAGAGGAAG CAGGATTCTTACAGAAAACAAGTGGT

This is the format which is understood by Swift's primer clipping tool, but in order to use other tools, I would like to convert it to BEDPE so I can then use standard tools to clip the primers from my bam files. However, I can't do that since I am not sure what the fields in this file are.
This is what I think thy are:

Chromosome
Target region start
Target region end
5' Primer name
5' Primer start
5' Primer end
3' Primer name
3' Primer start
3' Primer end
5' Primer sequence
3' primer sequence

Is my interpretation correct? Also, is this a standard file format or something specific to the Swift tool? I couldn't find any documentation about a "masterfile" format online.

Did you tried to contact with the tech support? Looking at the source code it seems that it accepts some "panel bedfiles"

@Llopis yes, I am talking to their tech support because I have other issues with the file. My question here is about this file format and whether it is some sort of standard format I am not aware of.

this is the swift standard primer file. You are a bit off in your interpretation of the columns (seems like you missed one as there are 12 columns, not 11 which threw everything off):
1) Chromosome
2) Target Start
3) Target Stop
4) Target Name (change which affected everything downstream in counting)
5) 5' Primer Start
6) 5' Primer Stop
7) 5' Primer Name
8) 3' Primer Start
9) 3' Primer Stop
10) 3' Primer Name
11) 5' Primer Sequence
12) 3' Primer Sequence
To convert to a BEDPE describing the primer locations, you only will need columns 1, 5, 6, 8, and 9 (and optionally 4). The last letter of the Region name for the denotes the strand orientation of the region. For example Region 1 has _F denoting forward or (+) strand for BEDPE usage, while region 2 has _R in name denoting reverse or (-) strand for BEDPE. Putting all this together, you can convert to BEDPE with:
$ awk -v OFS="\t" '{
            o=".";
            if($4~/_F\//){
                o="+"
            }
            else if($4~/_R\//){
                o="-"
            }
            print $1,$5,$6,$1,$8,$9,$4,".",o}' foo.masterfile
1   43814932    43814951    1   43815064    43815081    MPL_S505N-A519T_F/2 .   +
1   115256457   115256480   1   115256576   115256602   NRAS_G60-S65_new_R/1    .   -

Thanks, I meant 4 was the target name, I just wrote it wrong in the question. But why would I use the first 4 columns? I want the positions of the primers, not of the target. So, for a bedpe file, wouldn't I need columns 1,5,6,1,8,9 and then - or + depending on the "F" or "R"? Aren't the 1st 4 columns describing the target and not the primer?

I see, it was not clear to me that you wanted the location of primers and not the target. The primer orientation is also included in the primer names. The columns still remain the same, so if you would like primer info just use the information in the corresponding primer columns. The orientation of each primer is in its name as well. You would use columns 1, 5,6,7 for 5' primer and 1, 8,9,10 for 3' primer. Using name columns only for primer orientation (but BEDPE can also accept a 4th column for the name of the region if you'd like).

Yes, sorry, my question wasn't very clear on that (I edited now). Do you have any idea what the numbers at the end of the target names are? I mean the 2 at the end of MPL_S505N-A519T_F/2 and the 1 at the end of NRAS_G60-S65_new_R/1. Could they be the phase? And, if so, with respect to what?

By the way, I took the liberty of adding an awk command to make the conversion since I think that might be useful to future visitors. I hope you don't mind.

The awk command adds a nice touch! I am not super sure what the 1 and 2 mean with respect to the genes or orientations as with IDT orders I know they have a different meaning. Specifically, I think they refer to exon number but that could not be the case here. I don’t want to give you misinformation, but that is something you should check first with respect to your target set.


text:
Simulate and test CNV workflow?
I'd like to evaluate a CNV project.
My aim is to evaluate if the scripts are sufficient for calling reasonable CNVs. I know they have a paper, but their scripts may be buggy... and all papers are biased...
While there are CNV simulators on Google search, none seem popular in the community. How should I approach this problem?

Simulate CNV variants from a genome
Simple evaluation on the called variants against my simulation?

Did you try to follow the structure of their code? And using other metrics/programs to evaluate the CNV?

I cannot tell exactly from the post or the linked project what type of data these scripts are designed for, but I am going to assume that "Scripts involved in our workflow for detecting CNVs from WGS data using read depth-based methods" refers to short-read Illumina sequencing.
Simulations can be helpful, but they are usually imperfect. Nowadays we have long-read data available. I would suggest obtaining some long-read data for which short-read data is also available, doing CNV detection on both datasets independently, and comparing the results.
The problem with using long-read data is that any differences you find will be conflated with the biases of the technology. So, if you find a difference, then this approach will not tell you if the differences are due to the biases of the technologies, or due to bugs in the scripts. If you use simulations, then this conflation of biases remains.
Short of re-writing the pipeline, the only other suggestion I have is to borrow from software testing. Try constructing a few minimal examples that might also test some corner cases and see if the scripts behave exactly as expected.


text:
Reconstruct a protein's 3D structure from a pair-wise distance matrix
What is a good method to reconstruct a protein's 3D structure from a pair-wise distance matrix of residues?
I think the question is self-explanatory. I know of this method
but it seems that it uses contact maps, which actually use a certain cut-off. I don't have a cutoff based contact matrix but a matrix with actual distances.

You have distance between alpha carbons or beta carbons (or other atoms of the residues)? Usually a threshold is set because below certain measure it can't be and above certain threshold it doesn't mean a contact between the residues. Do you have any other information besides the distance matrix, like the sequence of the protein?

@llrs They are only carbon alpha distances and I want to reconstruct the backbone with carbon alphas and no more atoms.No more information, just the distances between carbon alphas for each residue.

You could use the multidimensional scaling approach, but I suspect it might propose invalid angles between the backbone. It would be possible to use it as a starting point before checking the angles between the backbone. But I don't know if there is any software to do this (it is a strange problem to have), or either you know the sequence or the sequence and the structure, but just the distance is strange.

Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.
But I don't enough about the context of your problem to tell if/how it is applicable.


text:
Does the ".full.aln" file produced by snippy-core contain all bases of my input sequences aligned to the reference genome?
I have a number of sequences and a reference genome.
I used snippy to align each individual sequence with the reference genome. I then used snippy-core on the *.aln files from snippy which produced a new *.aln file, and a *.full.aln.
It is my understanding that the new *.aln file contains a core SNP alignment with only bases shared by all of my sequences, and that the *.full.aln file contains the core genome alignment with all bases of each sequence aligned to the reference genome - regardless of variation between sequences. Is my understanding correct? If not where am I going wrong?
Further, the program snp-sites uses both the new *.aln and *.full.aln file to produce a file which can be run through FastTree to produce a phylogenetic tree. Does this just process the files to give the phylogenetic distance between each of my sequences? I don't really understand what this step achieves at all...

snippy-core is a tool that processes the snippy_outdir_{1..n} folders and produces the following files:

core.aln – a fasta format file that contains only the polymorphic sites (of {'A', 'C', 'T' and 'G'}) in your input set, whose length tends to decrease asymptotically as more isolates are added
core.full.aln – a fasta format file that contains all the core sites and base calls at any site where the reference has a site. The reference is not padded but the sequences mapped to it may be padded with gaps and 'N's to match the length of the reference.  The alignment length does not change as more isolates are added.

snippy{-core} was written by Torsten Seemann et al.  snp-sites was written by Andrew Page et al.   snp-sites will reduce any alignment to the core polymorphic sites (columns) and will still output a fasta.  Therefore the output of a core.aln processed through snp-sites should be the identical the input.  All fasta alignments can then be processed downstream using your tree inference package of choice (e.g., FastTree, IQ-Tree, MrBayes, BEAST ...).  If you want a distance matrix output from a fasta alignment, use something like snp-dists, also written by Torsten Seemann, or the APE package in R.  Be sure to consider how gaps and Ns are treated during this step.

Thank you so much for this. Was breaking my head for a long time. Its working now.


text:
bedtools single-nucleotide coverage in BED-specified regions for multiple BAMs
I have a BED6 (BED + name, score, strand info) file that defines some regions of interest.  I also have a set of BAMs corresponding to different samples.  I would like to obtain output similar to bedtools genomecov -d -5, which generates counts of read start positions, but I would like it restricted to regions of interest.
I can do this for individual BAMs and regions by filtering input with samtools view, e.g.,
bedtools genomecov -d -5 -strand '+' -ibam <(samtools view -b sample1.bam "chr1:12345-23456")

but this requires extracting position and strand info from the BED file, and then eventually I have to aggregate all the outputs back into a single file. This ends up being very slow (e.g., hours for ~50 regions of interest X 12 samples).
Is there a better way to accomplish this?
Ideally, something along the lines of how bedtools multicov works would be great:
bedtools multicov -bams /path/to/my/bams/*.bam -bed roi.bed

but it doesn't have the -d and -3/5 options to extract single-nucleotide read-end coverage.

Do you want coverage for each sample separately, or summed up across all samples?

50 ~1kb regions (as hinted at in the example) should not take hours to compute. Can you give more meta-information on the regions and BAMs that you are trying to process?

What's your end goal after making a text file with the values?


text:
Score predefined ChIP-seq peaks with MACS2 or equivalent
I have performed peak calling on a number of separate ChIP-seq experiments and would like to harmonise the peaks from each of the experiments in order to convert my data into a matrix for further analysis. I would like to do this by taking the union of the called peaks (available in BED / narrowPeak / MACS2 format) and retrieving the enrichment score / q-value for all of the experiments on each peak. Is this possible without writing my own extension to MACS2?

This is going to be an custom script so you're not going to touch MACS2 itself. The only question is how you want to summarize q-values from different samples. Other than that this is a job for bedtools/bedops.

Here's why I think I still need MACS - say Sample A has a peak at chrN:100-200, but Sample B does not. In order to avoid missing values, I would like to get the (not significant) enrichment score from that peak for Sample B using the MACS2 algorithm.

Ah, if you really want that then yes you'll need to use the (undocumented) MACS2 API.

You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like:
macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7

You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise.
But you can then use these "non-significant" values in combination with a bedtools map command, to map them to your predefined union of "actual" peaks.


text:
Fast and reliable alternatives to blast
After some unexpected results (and previously reported) I heard that there are other tools for finding similar sequences besides blast that are faster and more accurate.
I only found hmmer, but I don't know if it is faster and more accurate.
Is there any review comparing such tools? A quick google scholar search didn't report anything useful.
I want to search with short DNA sequences (~400), mainly interested for taxonomic imputation so I'll be looking for similar sequences and not distant homologs. I won't try to model genes and the sequence will be unspliced as it is from direct sequencing of a a region of the 16S gene.

Please [edit] your question and give us more details about what you will be looking for. Will these be nucleotide queries? Protein? Will you be looking for very similar sequences or distant homologs? Long sequences? Short ones? These details will affect your choice of tool.

My question is mainly theoretically, but I narrowed with the use case that prompted me to ask/think about this

The thing is you have tools like blat and megablat, FASTA (although I doubt anyone uses that anymore), hmmer, cd-hit, psi-blast, etc. It's really hard to answer unless you give a specific usage scenario. For instance, will you be trying to model genes? Would you want a tool like exonerate or genewise that can align a spliced cDNA sequence to a genome? Or will your 400 sequences be unspliced?

Ok, I realize this is more complicated than I thought. Thanks it is not my field/research focus :D. I updated the question with more information. Thanks for the guidance.

Try VSEARCH (or USEARCH). They are not BLAST replacement but can be a good choice for the right applications.

taxonomy assignment algorithms often use kmer-based methods which are faster than alignments (it looks like you are doing something like this?). You could just search for taxonomy assignment instead of alignment to find relevant alternative methods.

Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA

As I understand, the software tool Lambda is a viable, yet lesser known alternative to BLAST in the context of taxonomic classification of NGS data.

Why not shift the parameter inputs of Blast to exclude the species/subspecies which are being queried? Happy to show you how to do it, but that will pick up homologues. Alternative parse the blast output for e.g. 20000 outputs using E or identity within the threshold required (that's the way I'd do it).

All tools I’m aware of will be more-or-less as fast as blast. With kmer-based methods, the investment will be the up-front kmer database build.
Not sure if you’ve tried this, but NCBI has a 16S microbial blast database (lots of species, 16S rRNA only) just for this purpose. No point to search millions of nts for the 16S sequence every time if you know it’s SOME 16S sequence, so limit the search to only16S sequences.
I’d recommend you build a local blast database from ftp://ftp.ncbi.nih.gov/blast/db/16SMicrobial.tar.gz
BTW, another answer says blast won’t give you the alignment. Blast actually will give you alignments, you just have to ask it for that output format.

You can give a try to this tool: dnaservic.es
It's actually available experimentaly. It is a real local aligner, not just a search tool like BLAST. It really returns an alignment. It is usually slower than BLAST, but bot much slower in most cases. In contrast it usually returns alignments with better scores and not only partial alignments.


text:
What is a simple command line tool for doing Needleman-Wunsch pair-wise alignment on the command line
I have two DNA strings:

GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC

and

AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I want a tool that allows me to do something like this on the command line:
$ aligner  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

and receive an ASCII visualization of a pairwise alignment. Something like this would work:
GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---------
---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I prefer a tool that is packaged on bioconda. A tool that also does protein and RNA sequences is even more preferable.

Is this a toy example? If so, maybe update it to another one. These sequences are low complexity and the actual correct alignment output is unclear to me

The sequences are not important. Edits are welcome.

To get the alignment in your example, you have to use zero (or very small) gap cost at the ends of sequences. That is not the standard NW.

I'm not sure that's true? To quote from the EMBOSS package manual:
"In a Needleman-Wunsch global alignment, the entire length of each sequence is aligned. The sequences might be partially overlapping or one sequence might be aligned entirely internally to the other. There is no penalty for the hanging ends of the overlap. In bioinformatics, it is usually reasonable to assume that the sequences are incomplete and there should be no penalty for failing to align the missing bases."

You are looking for the needle program from the EMBOSS suite. Available in bioconda.
http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html
To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, extend_penalty 0.5) you'll need -auto. Thus your query above would be:
$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
#
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47

#---------------------------------------
#---------------------------------------

Install biopython with conda and use the following script:
#!/usr/bin/env python
import sys
from Bio import Align
aligner = Align.PairwiseAligner()
aligner.mode = "local"
alignments = aligner.align(sys.argv[1], sys.argv[2])
print(alignments[0])

The output is then:
.GGA-GGAGGGAG--AAGGAGGGAGGGAAGA-GGAGGGAG--AAGGAGGGAGGC
.|-|-||||||||--|||-|||-|||||-||-||||||||--|||-|||-|||.
AG-AAGGAGGGAGGGAAG-AGG-AGGGA-GAAGGAGGGAGGGAAG-AGG-AGG.

There are a LOT of alignments with that same score and you can tweak the details of the mode by changing aligner.mode.

yes, this should be pretty fast too since the alignment work is done in C (rather than python)

I'd just use something like t_coffee and write a little bash function to run it on the strings you give:
aligner(){
  tmp=$(mktemp);
  args=("$@")
  for ((i=0; i<$#; i++)); do
    printf '>%s\n%s\n' $i ${args[i]} >> "$tmp"
  done
  t_coffee "$tmp"  2>/dev/null | /bin/grep -Pv '^CLUSTAL|^\s*$'
  rm "$tmp"

}

If you add that to your shell's initialization file (e.g. ~/.bashrc if using bash on Linux, ~/.profile if using bash on macOS), you can run align2seqs with as many sequences as you want to align:
$ aligner 'GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG'
0               GG--AGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
1               nAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGG-GAAGAGGAGG
                    *********      *******  *********  * ** *  *
$ aligner 'GAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'GGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG' GAGGGAGGGAAGAGGAGGGAGAAGGATCCGAGGAAGAGGAGG AGAAGGAGGGAGGGAAGAGGCCCGAGATATATAAGGATCCGAGGAA
0               GA-GG----GAGAAGGAGGGAGGGAAGAGGAGGGAGAAGG--AGGGAGGC
1               GGAGG----GAGGGAAGAGGAGGGAGAAGGAGGGAG-GGA--AGAGGAGG
2               GA-GG----GAGGGAAGAGGAGGGAGAAGGATCCGA-GGA--AGAGGAGG
3               AGAAGGAGGGAGGGAAGAGGCCCGAGATATAT--AA-GGATCCGAG-GAA
                    *    ***      **   **     *       *    * *

Here's another option: https://github.com/noporpoise/seq-align; I've used it in the past as it is small and has no dependencies (just make and you're set).

This tool actually prompted my question! Do you know if anyone has packaged it?

Ah, no... I didn't notice the last part of your question.

Could you edit and include an example of how the OP could use this tool to do what they ask for? That will help future readers of this question.


text:
Finding the conserved region of protein in a given set of Kingdom
We have to find the conserved region of proteins in a given set of kingdom and we have to give more weightage to distantly related organism. What approach I should take to solve this problem?

Welcome to the site. In order to be able to answer the question better, could you explain what data/information do you have? Did you search for some terms related to your question? Why what you found didn't help? What else have you tried?


text:
Difference between genometools installed with anaconda vs apt-get?
I installed genometools with anaconda conda install -c bioconda genometools. Something is installed though apparently no binary called gt. If I type in gt Ubuntu suggests to install genometools via apt. What is the difference between both versions? Why is there no gt binary in the anaconda package?

These appear to be two separate collections of tools. The conda version is a collection of python scripts whereas the apt version is based on a C library. So, as you might expect, the python version does not have a gt binary.
They also have separate GitHub accounts conda genometools vs apt genometools.


text:
What is mate rescue in bwa mem?
BWA mem has the -S and -P tags for skip mate rescue and skip pairing; mate rescue performed unless -S also in use.
What do these do? I presume -P aligns read pairs independently of each other. Is that correct?
And what does -S do?

What version of BWA mem are you using?  Have you looked at the manual that gets installed with bwa?  "Please `man ./bwa.1' for the manual."

0.7.17; I did look at the manual page but the explanation is somewhat cryptic: "-P: In the paired-end mode, perform SW to rescue missing hits only but do not try to find hits that fit a proper pair."

Using -SP is equivalent to running bwa mem on each of the two mates as if they are single-end reads, but it formats the output as a proper paired-end output, i.e. with all pair-related flags added properly. Without -SP, by default bwa mem forces an alignment of a poorly aligned read if its mate is aligned somewhere. -SP turns off the forced alignment. We use -SP for mapping Hi-C reads.

If one read maps and the other doesn't, BWA attempts to rescue the other read by performing Smith-Waterman alignment with the unmapped mate, -S disables this mate rescue.
-P does indeed disables pairing, but mate rescues is still enabled - I think -P disables the setting of the proper pair flag. If you set -P and -S, BWA maps the paired reads essentially as single reads.

What does "attempts to rescue" mean? Why would the second read map using Smith-Waterman and not using the original algorithm?


text:
how to filter a multi-VCF for sex-specific genotypes?
I have a single VCF file of variants for a population that includes males & females and I would like to pull out all variants that are sex-specific. For example, sites where all males are 0/1 and all females are 0/0. Is there a way to do this using bcftools or vcftools? I could enumerate each individual name, or ask for all individuals named *.M and *.F.
I first tried running smoove on the male and female populations separately then comparing the two VCFs, but I'm working with a draft genome with loads of contigs. As only contigs containing variants are output in each VCF, and these are different between the two files, vcftools couldn't handle the comparison.
Someone here https://github.com/samtools/bcftools/issues/118 managed to get it working ike this:
$bcftools view -f PASS -i 'GT~"0/0.*0/1.*0/1.*0/1.*0/1"' combined.vcf.gz

but the same syntax gives me no variants (yet I can see some that meet my criteria in IGV).

Do you only want sites where either all males or all females are heterozygotes? Or do you also want sites that are more generally polymorphic in males or females but not the other?

Found a great tool with SnpSift!
From the examples: I want to keep samples where the genotype for the first sample is homozygous variant and the genotype for the second sample is reference:
cat variants.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] ) & isVariant( GEN[0] ) & isRef( GEN[1] )" > filtered.vcf

Given an input file input.vcf, I would take this approach:

Create two text files: one with the male sample names (males.txt) and one with the female sample names (females.txt) (one sample per line)
Create a BCF (for speed) for all the samples in each list while dropping the INFO field so that INFO/AC is recalculated (I'm not sure dropping the INFO field is strictly necessary). Then only keep sites with a minor allele count greater than zero:
bcftools view --samples-file males.txt input.vcf -Ou | bcftools annotate -x INFO  -Ou | bcftools view --min-ac 1:minor -Ob -o male_variants.bcf
bcftools view --samples-file females.txt input.vcf -Ou | bcftools annotate -x INFO -Ou | bcftools view --min-ac 1:minor -Ob -o female_variants.bcf

Use isec to create intersects and complements of the two variant files and store them in dir/*:
bcftools isec male_variants.bcf female_variants.bcf -p dir -Oz

Double-check dir/README.txt, but the files you seek should be named dir/0000.vcf.gz and dir/0001.vcf.gz
Disclaimer: This code is untested.

using vcffilterjdk  http://lindenb.github.io/jvarkit/VcfFilterJdk.html

or ask for all individuals named *.M and *.F.

java -jar dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().allMatch(G->(G.getSampleName().endsWith("M") && G.isHet()) || (G.getSampleName().endsWith("F") && G.isHomRef())); ' input.vcf

variant.getGenotypes().stream() get a stream of genotypes 'G' for the current variant.
allMatch all genotypes must match the statements below:
G.getSampleName().endsWith("M") && G.isHet() sample 's name ends with 'M' and is heterozygous
OR
G.getSampleName().endsWith("F") && G.isHet() sample 's name ends with 'F' and is Hom-ref

if you could provide seperate list of male and female samples names you can something like this:
$ bcftools view -e 'GT[@males.txt]!="0/1" && GT[@females.txt]!="0/0"' combined.vcf.gz

or if you don't want to take care of the exact genotype you can write it more generic:
$ bcftools view -e 'GT[@males.txt]!="alt" && GT[@females.txt]!="ref"' combined.vcf.gz

fin swimmer


text:
How to plot gene mutation frequencey in terms of percnatge alteration in samples?
I'm taking data from cbioportal of gene and their mutation frequency with percentage which I'm trying to categorize with various sub class.
The issue I'm facing is I have to plot it in mutation frequency which are already given. I'm not sure how to go about it.
This is my code (I'm attaching the data), but what I'm missing is the frequency count of patients. As an example NPM1 altered in 50 patient with 31.00% where the number of patient is 162, I'm going for stacked bar plot, as the frequency is already given, how can I just plot in terms of frequency ?
data
df.m1 <- reshape2::melt(ALL_NEW_RBP, id.vars = NULL)

library(dplyr)
d3 <- df.m1 %>%
  group_by(variable,value) %>%
  summarise(count=n()) %>%
  mutate(perc=count/sum(count))

colnames(d3) <- c("variable","value","count","percentage")

cbPalette <- c("black","darkblue","brown4","darkred","darkslategrey","deeppink","orange4","mediumorchid1","lightpink",
               "khaki1","gray70","cyan","coral1")

ggplot(d3, aes(x = factor(variable), y = percentage*100,
          fill = factor(d3$value,levels = c("3 ","4","5","6","7","8","9","10","11","12","13","31","NA","0")))) +
  scale_fill_manual(values = cbPalette, na.value = "lightsteelblue") +
  geom_bar(stat="identity", width = 0.7) +
  labs(x = "Categories", y = "percent alteration", fill = "Percentage") +
  theme_bw(base_size=15)+

  theme(axis.text.x=element_text(angle = 45, size=25, face="bold", hjust = 1),
        axis.text.y=element_text(angle=0, size=35, face="bold", vjust=0.5),
        plot.title = element_text(size=40, face="bold"),
        legend.title=element_blank(),
        legend.key.size=unit(1, "cm"),      #Sets overall area/size of the legend
        legend.text=element_text(size=30))

Could you attach the resulting image? Also, I think this question is more about how to use ggplot2 than about bioinformatics, you probably will have more answers if you post it on StackOverflow

added the figure..where the legend suggest the percentage of mutation you can have a look at the data and let know


text:
What i5 index should I use on the Illumina sample sheet for an unindexed p5 primer?
I have an upcoming run on a HiSeq X and most of the libraries in the pool have both i5 and i7 indices. However, some of the libraries were made with the IS4 p5 oligo and it is unindexed.
The IS4 indexing oligo is from the following paper. The sequence is below.

Meyer, Matthias, and Martin Kircher. "Illumina sequencing library
preparation for highly multiplexed target capture and sequencing."
Cold Spring Harbor Protocols 2010.6 (2010): pdb-prot5448.

                                 miseq index
                              5' TCTTTCC 3'
                                 >>>>>>>
 5' AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTT 3'
    ^------------------^         ^--------------------^
    flowcell binding seq        adapter-binding seq

When we sequence these libraries dual-indexed on the MiSeq, the machine reads into the adapter and returns the index sequence TCTTTCC.
Question
What will the index sequence be when the HiSeq reads it?
The more fundamental question is: How exactly do MiSeq and HiSeq machines read index sequences?
side note
Here is what an index p5 and unindexed p5 molecule look like:

Update: Sept 2019 - After running into this problem again, an Illumina representative agreed with Devon Ryan and said that setting up a run for single indexing is the best option when there are unindexed p5 ends on one or more libraries.

Illumina machines read the two indices in separate index sequencing reactions where they have sequencing barcodes specific for indices. If you have no index, then you should either not bother with creating either index read or mask that during demultiplexing. If you're mixing indexed samples with an unindexed sample, you can't really rely on any particular sequence being returned by the machine (i.e., just use the "undetermined indices" fastq file). The unindexed adapters tend to not have the full index sequencing primer binding site on them, so reads without index generally don't return much signal. In practice, they're likely to adopt barcode sequences of neighboring clusters. We see that a lot when we spike in PhiX. It has no barcode but you tend to get a small amount assigned to each sample either due to index hopping or (more likely) bleed through of signal from neighboring clusters (or small amounts of those clusters walking between nanowells and then their minority signal dominating during index sequencing).

Thanks for the suggestions! Now to see if the sequencing company will send me the bcl files. :)


text:
How to align output of grep --color=always? (To QC fasta/fastq files)
Grepping out short sequences from a fasta or fastq file is a really useful way to look at sequencing data. Using the option --color=always makes this even more useful, as you can visualize where the sequences appear in sequencing reads. For example:
zgrep 'GATCGATC' --color=always My_Hiseq_Data.fq.gz

To make this more useful, how can I align the output of grep --color=always such that the beginning of the highlighted parts are aligned?

Perhaps, grep is not the best tool to use in this case, but it should be in principle possible by using grep & sed. Here is an example showing three symbols around a match.
zcat My_Hiseq_Data.fq.gz | \
      grep -Eo '.{0,3}GATCGATC.*' | \
      sed -En 's/.*/   \0/; s/.*(.{3}GATCGATC.{0,3}).*/\1/p' | \
      grep --color=always GATCGATC

Here are some explanations:

zcat will just decompress a file.
grep will output all the matches with 3-symbol context (or less). This is needed because a line can have multiple matches.
first sed command will add 3 symbols before all the matches.
second sed command will output the sequence with 3 surrounding characters on the left and max 3 characters on the right. This is an equivalent of grep -Eo '.{3}GATCGATC.{0,3}'.
Finally, now that all matches are aligned, grep will color your sequence.

This code is possible to adapt to include more symbols of context.
I would recommend doing this in another language, though. E.g., in Python:
#!/usr/bin/env python3
import sys
import re

HL = '\033[93m'
ENDC = '\033[0m'

filename = sys.argv[1]
seq = sys.argv[2]
context = 10

seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (context, sys.argv[2], context))

for line in open(filename):
    line = line.rstrip()
    for m in seqre.finditer(line):
        print(m.group(1).rjust(context),
              HL,
              m.group(2),
              ENDC,
              m.group(3),
              sep='')

UPD
A more advanced Python version inspired by this answer with support for multiple files, gzip, etc.

#!/usr/bin/env python3
import sys
import re
import argparse
import gzip
from termcolor import  colored

def align_matches(f, seqre, context):
    # for every line
    for line in f:
        # get rid of newline character
        line = line.rstrip()
        for m in seqre.finditer(line):
            # for every match print right-aligned left-context,
            # match colored in red, and right-context
            print(m.group(1).rjust(context),
                  colored(m.group(2), 'red'),
                  m.group(3),
                  sep='')

# this function will check file extension to be able
# to open text files as well as gzip-compressed files
def myopen(fn):
    if fn.endswith('.gz'):
        return gzip.open(fn, 'rt')
    else:
        return open(fn, 'rt')

if __name__ == '__main__':
    # define script arguments
    parser = argparse.ArgumentParser(description='Aligned grep')
    parser.add_argument('pattern', type=str,
                        help='a valid regular expression')
    parser.add_argument('input', type=myopen, nargs='+',
                        help='files names, could be a gzip compressed file ending with .gz')
    parser.add_argument('--context', '-c', default=30, type=int,
                        help='number of context characters to show')
    args = parser.parse_args()

    # compile a regular expression with a match surrounded
    # by at most args.context number of characters
    seqre = re.compile(r'(.{0,%d})(%s)(.{0,%d})' % (
        args.context, args.pattern, args.context))

    # for every input file run align_matches
    for f in args.input:
        align_matches(f, seqre, args.context)

I wrote a Perl script to do this. Just like zgrep, it can search through as many files as you give it. You can use the -c switch to tell it how much context (how many characters) you want around the match. So -c10 will search for your target string and 10 characters on either side of it. If you don't use the -c, it will default to 155, which should include the entire sequence for most cases.
To use it, save as alignString.pl somewhere in your $PATH (~/bin should work if you're using Ubuntu, for example), make it executable with chmod a+x ~/bin/alignString and then run it like this:
alignString.pl GAAAATCAG file1.fastq.gz file2.fastq.gz file3.fastq.gz

For example:

The code:
#!/usr/bin/env perl

use warnings;
use strict;
use Term::ANSIColor;
use Getopt::Std;

my %opts;
getopts('c:',\%opts);
my $context = $opts{'c'} || 50;

## This subroutine does the alignment
sub align{
  ## Read the line and the target pattern
  my ($line, $pat) = @_;
  ## This will find 0 to $context characters on either side
  ## of the search pattern. Note that it will find the first
  ## occurrence of the pattern.
  $line =~ /(.{0,$context})$pat(.{0,$context})/;
  ## Save the strings to the left and right of the pattern
  my $left = $1;
  my $right = $2;
  ## If the left hand string is shorter than the value of $context,
  ## pad it with spaces so it can be aligned.
  while (length($left) < $context) {
    $left = " " . $left;
  }
  ## Return the padded left hand string, the pattenr colored red
  ## and the right hand string unchanged.
  return "$left" . colored($pat,"red") . "$right";
}

## The search pattern is the 1st argument
my $pattern = shift;
## Any remaining arguments are taken as target files to search
my @fastqFiles = @ARGV;

## Run the zgrep command
open(my $grepStream, '-|', "zgrep $pattern @fastqFiles");
## For each line returned by zgrep
while (my $line = <$grepStream>) {
  ## Align this line
  $line = align($line, $pattern);
  print "$line\n";
}

wonderful! I've resisted learning perl for a long time but it seems useful sometimes...

Beatiful! :) I didn't know about Term::ANSIColor, it is so simple to use! Your left padding while loop can probably be shortened to: $left = " " x ($context - length($left)) if length($left) < $context;

@PeterMenzel thanks :) and yes, various bits of this can be golfed but I was trying to keep it as clear to non Perl people as possible.

@terdon - over a year later I found an edge case where this doesn't work. If you search for a homopolymer like 'TTTTT', if there is a case where there are multiple matches (such as in a longer homopolymer 'TTTTTTTTTT'), this perl script will align to the rightmost match- not the leftmost


text:
Fast filtering of intervals not falling within a certain distance from known genes
I would like to filter a bed file with intervals, ie. in the format of:
chr1    13800   14301
chr1    15500   16001
chr1    19400   19901
chr1    22800   23301

In particular, I want to filter out intervals that fall far from known genes, let's say > 1KB.
I was thinking of querying each interval (in parallel) and check through a known set of coordinates from a gtf, e.g. from Homo_sapiens.GRCh38.93.gtf by filtering for gene coordinates:
1       havana  gene    168317497       168376876       .       +       .       gene_id "ENSG00000227722"; gene_version "2"; gene_name "AL022100.1"; gene_source "havana"; gene_biotype "transcribed_processed_pseudogene";
1       mirbase gene    168375524       168375621       .       +       .       gene_id "ENSG00000207974"; gene_version "1"; gene_name "MIR557"; gene_source "mirbase"; gene_biotype "miRNA";
1       havana  gene    168400829       168495685       .       -       .       gene_id "ENSG00000228697"; gene_version "2"; gene_name "AL023755.1"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168401483       168407274       .       -       .       gene_id "ENSG00000283255"; gene_version "1"; gene_name "AL022100.2"; gene_source "havana"; gene_biotype "lincRNA";
1       havana  gene    168449672       168451211       .       -       .       gene_id "ENSG00000215833"; gene_version "3"; gene_name "QRSL1P1"; gene_source "havana"; gene_biotype "processed_pseudogene";

However, since the bedfile is quite big I am looking for the fastest solution since checking each interval may be quite expensive.

Here's a way to use BEDOPS, which was designed to work fast by using sorted input. Other tools now use sorting to accomplish similar performance benefits.

Convert GTF annotations to a sorted BED file of genes:
$ awk '($3=="gene")' annotations.gtf | gtf2bed - > genes.bed

Sort your intervals, if unsorted:
$ sort-bed intervals.unsorted.bed > intervals.bed

Filter intervals against 1kb-padded genes:
$ bedops --merge genes.bed | bedops --range 1000 --everything - | bedops --not-element-of 1 intervals.bed - > intervals.filtered.bed

We first --merge genes to simplify dealing with overlapping genes and further speed downstream filtering.
The file intervals.filtered.bed will contain intervals from intervals.bed which lie outside of a 1kb-window around merged gene space.

To run all of this as a one-liner (assuming use of the bash shell):
$ bedops --merge <(awk '($3=="gene")' annotations.gtf | gtf2bed -) | bedops --range 1000 --everything - | bedops --not-element-of 1 <(sort-bed intervals.unsorted.bed) - > intervals.filtered.bed

This creates the same output result, but avoids creation of intermediate files and so should run even faster.

Note: BEDOPS is agnostic about chromosome sizes — this toolkit works with any genome assembly. However, using --range takes care of ranges that would have a negative start position, setting their start coordinates to 0. Bounds aren't really relevant here, in any case, because the bedops --not-element-of step queries input intervals, which should already be within chromosome bounds, unless they were generated incorrectly.

I think there's a missing quote on the first awk call

When executing 1 I get Error: Potentially missing gene or transcript ID from GTF attributes (malformed GTF at line [1]?)

https://www.biostars.org/p/206342/#206350

In one line, using bedtools
zcat  Homo_sapiens.GRCh38.93.gtf.gz \
 | awk '$3=="gene"' \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

This first takes the genes and filters them for the third column being gene. This is important because all GTF lines contain the word gene as gene_id is a manditory attribute for GTF.
Next it passes it to slop, which adds 1 kb either side (ensuring not to go below 0 or off the end of the chr: this isn't that important because we are never outputting these regions anyway, but bedtools requires it). We tell bedtools to work on the standard in by using -i -.
Finally we intersect these gene regions with the intervals.bed. -u says output the line from -a if it overlaps with at least one entry in -b and passing -b - means read file b from the standard in.
This will only work if your GTF contains gene lines, which the GENCODE/Ensembl ones do. However, it is interesting to note that these GTFs are not actaully valid GTF files, as to be valid a GTF file must have both a gene_id and transcript_id attribute, and the gene lines in the GENCODE GTFs to not have this.
If your GTF does not contain gene lines, but can create them by merging all the exons in all the transcripts and filling in the gaps using gtf2gtf from cgat.
 cgat gtf2gtf --method=merge-transcripts --stdin Homo_sapiens.GRCh38.93.gtf.gz \
 | bedtools slop -b 10000 -g contigs.tsv -i - \
 | bedtools intersect -u -a intervals.bed -b -

First I prepared a bed file in which the gene intervals are augmented by 1KB before and after the gene start and end coordinates. Then I intersected this bed file with my original one with the option -wa, therefore retaining only the intervals in my original bed file that intersect with the bedfile produced from the original gtf, filtered by gene regions. The -u flag is also used to print only one line per overlap between the two files, in case an interval in the first overlaps with multiple intervals in the second.
# download information about chromosomes
wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.chrom.sizes ./

# Original GTF has a different notation than UCSC
cat Homo_sapiens.GRCh38.93.gtf | \
  sed -e $'s/^./chr&/' -e $'s/^chrMT/chrM/' | \
  awk '$3=="gene"' | \
  bedtools slop -b 1000 -g hg38.chrom.sizes -i - | \
  sed -e $'s/\tgene_id.*gene_name/\t/' -e $'s/gene_source.*//' | \
  tr -d '"; ' | \
  awk -v OFS=$'\t' '{print $1, $4, $5, $9}' > genes_1kb.bed

bedtools intersect -wa -u -a intervals.bed -b genes_1kb.bed > intervals.1kb.bed

note that this will filter out intervals that overlap a gene, or 1kp upstream, but will leave features that are within 1kbp downstream. And it will only work if the GTF has gene lines, which many GTF files don't have.

@IanSudbery, why this will filter out intervals that overlap a gne or 1kb upstream?

Oops I noticed now I forgot to add 1kb to the end. Thanks! Does it solve the problem?

@gc5 your awk command could have been replaced with bedtools slop ref which has the added advantage of making sure the intervals do not extend more than the chromosome length. It also reduces the chances of zero or one base coordinate errors (tried and test tool).

@fridaymeetssunday amazing! I didn't know it, thanks


text:
API/website for blasting peptide sequence against database of all plant proteomes
I would like to blast a peptide sequence against all available plant proteome databases, i.e. a blastp metaproteome analysis.
What tools are available? Which database can I use?

Welcome to the site lars. You seem aware of blastp, did it not solve your problem? Or do you expect something else of blastp? Could you explain what have you tried?  (to avoid suggesting the same you have already done)

Thank you for the very quick reply. So far I installed ncbi-blast-2.7.1+. With makeblastdb I generated *.phr, *.pin, *psq from a fasta. I now try blastp against this database. Next steps would be to compile a metadatabase from all plant proteome databases I can find.
I wonder if all of this work has already been done. Is there a website/API that I can use? I simply want to find all plant proteins which contain the amino acid sequence TEKKRRETVEREKE. Thanks.

Please [edit] the question to include this information. In addition, what sequences are in your fasta file? I think that you can download fasta files from NCBI. You  could download the sequences from plants proteins and make the database from that (if this is not already what you did)

At the moment, I work only with arabidopsis thaliana. Is there a plant meta-database? My main objective is to figure out in which plant species this peptide is expressed.

For example, PPDB is only for arabidopsis thaliana and maize.

Sorry I don't work with plants (nor with blast) and I can't answer where to find such data. But add this relevant information to the question, so make it easier to those who might be able to answer you.

A simple BLAST on the UniProt website does the trick (target database plants).
Best matches are subsequences ETVEREKE and EKKRRETV in Green alga and Eelgrass respectively.


text:
cDNA and alignment mapping
I am confused with RNA seq alignment. My understanding is that after Mature mRNA is isolated from the cell, it is then fragmented and using reverse transcriptase enzyme a cDNA copy is created which is further amplified. This cDNA is then mapped with the reference genome. If the cDNA does not have introns how does it map with the reference genome.
Can someone please help

tl;dr
Spliced aligners just make many small matches along the read pair or single read to determine splice junctions and output an alignment for the whole read, as long as it matches somewhere in the genome.
mRNA and translation
In the cell, after transcription, introns are spliced out of RNAs to become mature RNA. The protein-coding sequence does have sequence that is chimeric from cis-located portions of the same chromosome that do not actually exist.
              |<<exon 1>>|  intron    |<<exon 2>|
Genome: ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

mRNA  : ----ABCDEFGHIJKLYZ123456789aaaaaaaaaaa

As you said, if we only even consider small stretches of DNA, the sequence JKLYZ does not exist in the genome, only in the mRNA. So, we need a special way to map mRNA-seq reads to the genome that correctly splits at splice junctions.
How do we map mRNA-seq reads to a genome?
As a result, parts of the mRNA-seq reads will align well to one part of the genome, but only for part of the read. This requires a different alignment algorithm that considers multiple local high-score alignments for different parts of the read.
For example, see how a genomic aligner versus a splice aligner might treat the same mRNA-seq read on our example genome.
                         |<<exon 1>>|  intron    |<<exon 2>|
Genome:            ------ABCDEFGHIJKLmnopqrstuvwxYZ123456789------------------

spliced aligner  : ------ABCDEFGHIJKL------------YZ123456789
genomic aligner  : ------ABCDEFGHIJKL

There is a possibility that the genomic aligner will not recognize that the second half of the read is an exon due to the large intron gap, and it may clip the alignment. Spliced aligners just take this into account and output several alignments for the same read regardless of gap size. The STAR algorithm, for example, uses many small matching pieces of a single read to determine where the splice junctions are. HiSat2 uses small approximate matches using a FM-index to find seeds for spliced alignment.


text:
peakranger bcp segmentation fault
I'm running into a segmentation fault running peakranger's bcp on Arch Linux. Binaries were compiled from source just today. Any idea how I can prevent this?
$ peakranger bcp -d chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam -c chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam --format bam -o chip_seq/align_count/Sample_DNA_loA/call_peaks/bcp --gene_annot_file /data/sai/annots/musmus_GRCm38.94.gtf.gz --verbose
program version:          1.18
Data files:
 File format:             bam
 Sample file:             chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam
 Control file:            chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam
Qualities:
 P value cut off:         0.0001
 sliding window size:     500
 Read extension length:   200
Output:
 Regions:                 chip_seq/align_count/Sample_DNA_loA/call_peaks/macs2_region.bed
 HTML reports:            Enabled
 Plot region length:      6000
 Annotation file:         /data/sai/annots/musmus_GRCm38.94.gtf.gz

Reads statistics:

 Treatment reads +:       2125600
 Treatment reads -:       2113802
 Average read length:     20
 Control reads +:         7582218
 Control reads -:         7575995
 Average read length:     20
 Verifying reads...
Warning: JH584295.1 only contains reads in the positive strand. The chromosome is removed.
Warning: GL456382.1 only contains reads in the negative strand. The chromosome is removed.
Warning: No reads were found in GL456382.1 of the treatment dataset. The chromosome is removed.
Warning: No reads were found in JH584295.1 of the treatment dataset. The chromosome is removed.

Reads statistics after correction:

 Treatment reads +:       2125597
 Treatment reads -:       2113800
 Control reads +:         7582196
 Control reads -:         7575947

 Calling peaks...

Segmentation fault (core dumped)

Reads were trimmed with trim galore and aligned with hisat2.
Update 2018-10-29: additional information of input and machine
$ samtools flagstat chip_seq/align_count/Sample_DNA_loA_h3k4me3/alignment/Sample_DNA_loA_h3k4me3.hisat2.bam
33521569 + 0 in total (QC-passed reads + QC-failed reads)
946323 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
5185725 + 0 mapped (15.47% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)
$ samtools flagstat chip_seq/align_count/Sample_DNA_loA_input/alignment/Sample_DNA_loA_input.hisat2.bam
22508880 + 0 in total (QC-passed reads + QC-failed reads)
3359461 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
18517674 + 0 mapped (82.27% : N/A)
0 + 0 paired in sequencing
0 + 0 read1
0 + 0 read2
0 + 0 properly paired (N/A : N/A)
0 + 0 with itself and mate mapped
0 + 0 singletons (N/A : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

I tried this with two other input bam files and also hit a segfault.
Machine is a 36-core (72 hyperthread) workstation with 500GB of RAM.
$ peakranger --version
1.18
$ g++ --version
g++ (GCC) 8.2.0
$ uname -r
4.18.3-arch1-1-ARCH
$ pacman -Q boost
boost 1.67.0-7
$ pacman -Q zlib
zlib 1:1.2.11-3

Thanks @Llopis, I've clarified the question and added information on the bam files and the machine.