Contigs with no read coverage (zero reads mapped) after mapping using Bowtie2

I recently mapped reads using Bowtie2 against viral genome segments to find the read counts for each of the segments. To my surprise, I found contigs with no reads mapped on them.

That is when I came to know that bowtie2 used END to END read alignment using default settings (see here). In reality, the contigs from many assemblers are further broken down into kmers and are stitched together. In this case, I used the Megahit metagenome assembler, which uses the default sequence of 21, 41, 61, 81 and 99 kmers (Reference). So, aligning reads using the END to END parameter settings may not work sometimes resulting in zero reads count or no read coverage for the particular contig. 

In the figure down, we see the PB2 and NP (colored read) using END to END has no reads mapped but the contig is created by the assembler. This more likely seems to be happened because the assemblers use kmer approach to create contigs. This is when, changing the setting from END to END to LOCAL makes more sense. 

In the local alignment when some of the bases at the ends of the read do not participate, they are omitted (or "soft trimmed" or "soft clipped") from the beginning or from the end. That's how we see PB2 and NP have now read counts. For other segments, the read counts seemed to increase. 

So, in this particular case of alignment, local alignment of the reads using bowtie2 makes more sense.

For more discussion, see here:

• What Could Be The Situation That Contigs With Zero Reads And Singleton Contigs Are Included In The Assembly?
• Bowtie2 parameters for best alignment
• Bowtie2 and BWA-MEM giving very different results in metagenomic data

Downloading viral genome database for viral read classification - Metagenomics

1.  Kraken2 Metagenomic Virus Database - redirects to globus - does not have an option to download
   
   
2. Default kraken2 command:
• kraken2-build --download-library viral --db $DBNAME
  
  Results in the error - adding --use-ftp did not help
  rsync: getaddrinfo: ftp.ncbi.nlm.nih.gov 873: Name or service not known
  rsync error: error in socket IO (code 10) at clientserver.c(127) [Receiver=3.1.3]
  Error downloading assembly summary file for viral, exiting.
  
  
• Tried changing code in specific scripts based on this error thread - still no luck! :(
• Tried changing code in specific scripts based on this error thread - still no luck! :(
  
  
3. A python script that helps with updating the kraken databases :
• error - FileNotFoundError: [Errno 2] No such file or directory: 'assembly_summary_refseq.txt'
  
  
4. Downloaded NCBI Viral database from here - have not tested - but fasta sqeuences needed to converted to kraken2 database format.
   
   
5. Downloaded Viral RefSeq database from a PeerJ paper - worked!
• Pre-compiled databases
• Looks comprehensive! Size is big (6.6 gb)
• Could classify the viral reads using kraken2 command

Resolved-How to load an HTML page on Github as a normal HTML page we see in browser?

Github is all about code and there is no automatic way to load and view a .html file uploaded in our github repo. So, only the html code is displayed when clicked on the .html document.

A simple way to overcome this problem is to prepend the following before the actual link.

For example, 

my actual html file is located at:

https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.html

to load it as html file, we should prepend :

https://htmlpreview.github.io/?

Finally, my link in the browser should look like this:

https://htmlpreview.github.io/?https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.html

BRIG genome visualization missing rings - Solved!

Problem: 

Missing rings for the query sequences and BRIG doesn't give any error message.

Solution that worked for me:

This is an issue with java version. I was running Java 1.8 version and had same problem as the OP's question. Then I followed the answers and finally changing the java version worked for me.

I ran the following steps

1) Go to link: https://www.oracle.com/java/technologies/javase-java-archive-javase6-downloads.html 

2) mkdir Java1.6

3) Download "jdk-6u45-linux-x64.bin" and save in Java1.6

3) cd Java1.6 && chmod +x jdk-6u45-linux-x64.bin

4) You should see java file in Java1.6/jdk1.6.0_45/bin

5) Java1.6/jdk1.6.0_45/bin/java -version - should give "1.6.0_45"

6) Open BRIG using Java1.6 as follows:

Java1.6/jdk1.6.0_45/bin/java -jar BRIG.jar

Note: I did not uninstall my original Java1.8 rather used the downloded java1.6 execulatable file using path. That way I retain both versions without any problem.

This gave me all the rings without any issue!! Hope this helps someone!

https://www.researchgate.net/post/Why_does_BRIG_fail_to_create_rings_for_some_of_the_genomes

Posted my answer in research gate - click above link to see the discussion

SRST2 requires a specific naming convention for input fastq files

 SRST2 Trouble installing in server

 Also wanted to mention that SRST2 has been there for a while now (last updated was 2015) and I had some issues with running server because of dependency issues of bowtie2 tool (libstd, libtbb). 

 

SRST2 Trouble with specific versions of bowtie

 

Also, SRST2 only needs a specific version of bowtie2 (between 2.10 to 2.2.9 - but the latest is 2.4.1) causing some problem to work on server. 

 

SRST2 installation

 

Finally installed SRST2 on using conda. https://anaconda.org/bioconda/srst2

 

SRST2 requires a specific naming convention

 

An other issue with SRST2 is the naming, since it accepts the names like _S1_L001_R1_001.fastq.gz. 

For those samples which does not have a convention " _S1_L001_R1_001.fastq", created a symlink dummy with " _S1_L001_R1_001.fastq". 

For example, "filename_L001_R1_001.fastq" does not follow the SRST2 input convention, so created a symlink file filename_S1_L001_R1_001.fastq.

 

Nextflow: Missing 'bind' declaration in input parameter - resolved

$ time nextflow test4_20200322.nf

N E X T F L O W  ~  version 21.04.3

Launching `test4_20200322.nf` [happy_leakey] - revision: d7d5d88709

reads: /data01/nextflow_test/2022_test/sample_data/*_{1,2}.fastq

Missing 'bind' declaration in input parameter

 -- Check script 'test4_20200322.nf' at line: 9 or see '.nextflow.log' file for more details

Paste this line at the 2nd line : nextflow.enable.dsl=2

Tips for uploading data to NCBI SRA

 ### Case: 

So, If have 20 bacterial samples, then I will have 20 x 2 = 40 fastq files for Illumina paired data and 20 x 1 = 20 fastq files for Nanopore data that I need to upload.

1) Create a Bioproject (Get the Bioproject accession)

2) Create a Biosample (Use Bioproject accession to fill the Microbe template (template depends on the organism bacterial in this case;can be downloaded from NCBI only as going through the steps)

3) For Illumina data 

   - Create a SRA submission 1

   - Need to fill the metadata template using the SAMN numbers generated from Biosample submission.

   - Upload the data either using aspera or FTP.

4) For Nanopore data 

   - Create a SRA submission 2

   - Need to fill the metadata template using the SAMN numbers generated from Biosample submission.

   - Upload the data either using aspera or FTP.

# Pointers:

1. Always best to create a bioproject first and get PRJNA number

2. Then create a BioSample project - in the Microbe1.0 excel provide the Bioproject (PRJNA) number

3. Then do the SRA submission.

4. For the SRA submission, if you have two types of data such as Ilumina and Nanopore sequences for the same sample, then in the SRA_meta_acc file have the Illumina information and followed by nanopore details. 

5. The SAM number (Biosample number) will be same for both Illumina and Nanopore data

6. If your data samples are more than 1000, then single submission is not possible in NCBI. So you have to split the data accordingly.

7. aspera command line seems a bit confusing to use, ftp seems to be straightforward. Just need to get username password from SRA submission

8. Advantage with ftp is: Even if sometimes the internet connection drops, we can resume from where transfer stopped

9. We cannot upload half of our data with FTP and half of our data aspera. It is either aspera or FTP. 

Extract accessory genes as fasta files from the roary output for downstream functional characterization

1. First generate the locus_tag of the accessory genes (output file: pan_genome_results)

query_pan_genome -a complement -g clustered_proteins *.gff
mv pan_genome_results accessory_genome_locus_tags.list

2. Use the locus_tags of the accessory genes to extract the sequences from prokka folder

fgrep 'locus_tag=' *.gff >Combined_local_gffs.txt # search 'locus_tag' in all gff file and combine - for speeding purpose

time for locus_tag in $(awk '{print $2}' accessory_genome_locus_tags.list); do 
file=`fgrep -m1 "$locus_tag" Combined_local_gffs.txt | awk -F':' '{print $1}' | sed 's/.gff//g'`; # extract the samplename (foldername) where the gene is present
awk -v var="$locus_tag" 'BEGIN {RS=">"} $0~var {print ">"$0}' Prokka/MySampleFolders/"$file"_prokka_out/"$file".ffn >accessory_gene_representatives/"$locus_tag".fasta;  
done 

These genes can be combined and can be functionally characterised for downstream analysis

Extracting 5477 genes took me around ~19 mins by this way

Add file name after ">" in a multi fasta file

More solutions here: https://www.biostars.org/p/435813/

awk '/>/{sub(">","&"FILENAME"_");sub(/\.fasta/,x)}1'  assembly.fasta | sed 's/ /_/g' | grep '>'

Dowloading raw fastq data from European Nucleotide Archive (ENA) using ERR IDs

 Step1: Paste the comma separated ERR IDs in the address bar like this:

https://www.ebi.ac.uk/ena/browser/view/ERR3417384,ERR3418576,ERR3418577,ERR3307235?show=reads

Step 2: This will list the ID in the webpage like this:

Step 3: Click Download XML file - this downloaded ena_data_20211215-0958.xml file on to my desktop

Step 4: In the linux terminal:

$ fgrep 'filename="run' ena_data_20211215-0958.xml | 

sed -e 's/.*run/ftp.sra.ebi.ac.uk\/vol1\/run/g' -e 's/\".*//g' >downloadlinks.list

Step 5 : Looping through the links and using wget command from linux, we can download like this:

$ time for d in $(cat downloadlinks.list); do echo $d; wget $d; done

Step 6 : But, wget command was very slow and reconnecting many times like this:

--2021-12-16 09:50:41--  (try:13)  http://ftp.sra.ebi.ac.uk/vol1/run/ERR329/ERR3297045/1_P_PA_1.fastq.gz
Connecting to ftp.sra.ebi.ac.uk (ftp.sra.ebi.ac.uk)|193.62.197.74|:80... connected.
HTTP request sent, awaiting response... 206 Partial Content
Length: 1842088178 (1.7G), 606664946 (579M) remaining [application/octet-stream]
Saving to: ‘1_P_PA_1.fastq.gz’

1_P_PA_1.fastq.gz                                     67%[+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++                                       ]   1.15G  --.-KB/s    in 15m 0s  

2021-12-16 10:05:43 (0.00 B/s) - Read error at byte 1235423232/1842088178 (Connection timed out). Retrying.

--2021-12-16 10:05:53--  (try:14)  http://ftp.sra.ebi.ac.uk/vol1/run/ERR329/ERR3297045/1_P_PA_1.fastq.gz
Connecting to ftp.sra.ebi.ac.uk (ftp.sra.ebi.ac.uk)|193.62.197.74|:80... connected.
HTTP request sent, awaiting response... 206 Partial Content
Length: 1842088178 (1.7G), 606664946 (579M) remaining [application/octet-stream]
Saving to: ‘1_P_PA_1.fastq.gz’

1_P_PA_1.fastq.gz                                     67%[+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++                                       ]   1.15G  --.-KB/s    in 15m 0s  

2021-12-16 10:20:54 (0.00 B/s) - Read error at byte 1235423232/1842088178 (Connection timed out). Retrying.

--2021-12-16 10:21:04--  (try:15)  http://ftp.sra.ebi.ac.uk/vol1/run/ERR329/ERR3297045/1_P_PA_1.fastq.gz
Connecting to ftp.sra.ebi.ac.uk (ftp.sra.ebi.ac.uk)|193.62.197.74|:80... connected.
HTTP request sent, awaiting response... 206 Partial Content
Length: 1842088178 (1.7G), 606664946 (579M) remaining [application/octet-stream]
Saving to: ‘1_P_PA_1.fastq.gz’

1_P_PA_1.fastq.gz                                     81%[+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++================>                      ]   1.40G  --.-KB/s    in 26m 54s 

So, used aspera pluggin to download at better speed.

Step 7:  Locate asperaweb_id_dsa.openssh in ubuntu. The absolute path is required for the next step

$ locate asperaweb_id_dsa.openssh
/home/user/.aspera/connect/etc/asperaweb_id_dsa.openssh

Step 8: parsing the downloadlinks.list and creating a new file for aspera download.

$ cat downloadlinks.list | sed 's/ftp.sra.ebi.ac.uk\/vol1\/run\///g' >downloadlinks_ascp.list

$ time for d in $(cat downloadlinks_ascp.list); do echo $d; ascp -v -l 300m -P33001 -i /home/prakki/.aspera/connect/etc/asperaweb_id_dsa.openssh era-fasp@fasp.sra.ebi.ac.uk:vol1/run/$d /ena/download/folder/ ; done