Installing pbh5tools on OSX yosemite

I need to assemble PacBio amplicons of olfactory receptors. One of the main challenges I have been having is that most of the PacBio tools (SMRT-analysis, Celara assembler) focus on assembling whole genomes. Using an approach from Larson, et al. (2014) in BMC Genomics, I am going to do quality filtering using pbh5tools that will filter the reads at a minimum of passes (e.g. 4; though my mean number of passes is quite high for the data (~20)).

To install pbh5tools, I dealt with a series of challenges on Yosemite:
Download from github.

cd pbh5tools-master/

sudo python setup.py install

I ran into this pesky error.

In file included from /System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/core/include/numpy/ndarraytypes.h:1760:

/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/numpy/core/include/numpy/npy_1_7_deprecated_api.h:15:2: warning: 

      "Using deprecated NumPy API, disable it by "          "#defining

      NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION" [-W#warnings]

#warning "Using deprecated NumPy API, disable it by " \

 ^

In file included from /tmp/easy_install-_W_eCV/h5py-2.5.0/h5py/defs.c:279:

/tmp/easy_install-_W_eCV/h5py-2.5.0/h5py/api_compat.h:27:10: fatal error: 

      'hdf5.h' file not found

#include "hdf5.h"

         ^

1 warning and 1 error generated.

error: Setup script exited with error: command 'cc' failed with exit status 1

First I tries to install hdf5 tool packages and then export the path there, but that wasn't working.

From the setup.py, I see that it needs:
    install_requires=[
        'pbcore >= 0.8.0',
        'numpy >= 1.6.0',
        'h5py >= 1.3.0'
        ]

sudo pip install pbcore

sudo pip install numpy

sudo pip install h5py

Then,

sudo python setup.py install

...

Using /Library/Python/2.7/site-packages

Searching for pysam==0.8.3

Best match: pysam 0.8.3

pysam 0.8.3 is already the active version in easy-install.pth

Using /Library/Python/2.7/site-packages

Finished processing dependencies for pbh5tools==0.8.0

Voila!

Is deduplication for you?
Re: @brianbushnell on seqanswers:Deduplication of reads is possible with dedupe, but it requires a lot of memory (~1kb per read) so it is only practical with HiSeq data if you have high-memory nodes. If you have a reference, deduplication based on mapping will use much less memory, though it will also be much slower. And whether or not deduplication is a good idea depends on the experiment; i.e., probably not for quantitative RNA-seq, and also probably not for unamplified libraries in general.

Suggested order:
1. Initial quality control BGI did this
2. Quality and/or length-based reads discarding; trimming/discarding of N-containing reads BGI did this
3. Adapter removal BGI did this
4. phix/contaminant removal
5. [optional] Error-correction --use ecc.sh
6. [optional] Deduplication (depending on experiment) --use dedupe.sh
7. Merging of PE reads --use bbmerge.sh
8. Quality trimming + phix/contaminant removal --use bbduk.sh
9. Contamination check

More advice:
As for the order of quality trimming and merge, they could be done either way. But, BBMerge internally performs quality-trimming for reads that don't merge prior to quality trimming, so for that particular program, quality-trimming after merging tends to be better.

Note that after merging, you will have 2 sets of reads (merged and unmerged) that must be processed independently, since not all of them will successfully merge

Also useful reading about merging your paired end reads before assembly:
http://thegenomefactory.blogspot.com/2012/11/tools-to-merge-overlapping-paired-end.html

Okay, time to try...

Step 1: Quality trimming

#quality control trimming of paired-end reads

#DR086_E.bombifrons_MOE quality control and plots

bbduk.sh -Xmx1g in1=DR086_Erbom_MOE_1.fq in2=DR086_Erbom_MOE_2.fq out1=DR086_Erbom_MOE_1_clean.fq out2=DR086_Erbom_MOE_2_clean.fq minlen=25 qtrim=rl trimq=10 qhist=qc_plots/DR086_Erbom_MOE_qhist.txt aqhist=qc_plots/DR086_Erbom_MOE_aqhist.txt lhist=qc_plots/DR086_Erbom_MOE_readLength.txt bqhist=qc_plots/DR086_Erbom_MOE_bqhist.txt overwrite=true

About 8-9% of each read was trimmed.

Normal level sequence duplication in mammals is 20 million reads

Normal sequence bias at beginning of reads due to nonrandom hybridization of random primers.

#I put this all in a shell script

sh BGI_qc.sh >log.txt

I tried to do error correction but these files are HUGE!! Need to remove duplicates first.

sh ecc.sh in1=DR086_Erbom_MOE_1.fq in2=DR086_Erbom_MOE_2.fq out1=DR086_Erbom_MOE_1_corrected.fq out2=DR086_Erbom_MOE_2_corrected.fq

/usr/local/bin/bbmap//calcmem.sh: line 111: [: unlimited: integer expression expected

java -ea -Xmx-410m -Xms-410m -cp /usr/local/bin/bbmap/current/ jgi.KmerNormalize bits=16 ecc=t passes=1 keepall dr=f prefilter in1=DR086_Erbom_MOE_1.fq in2=DR086_Erbom_MOE_2.fq out1=DR086_Erbom_MOE_1_corrected.fq out2=DR086_Erbom_MOE_2_corrected.fq

Invalid maximum heap size: -Xmx-410m

Could not create the Java virtual machine.

#even when I increased the memory, no cigar

sh ecc.sh in1=DR086_Erbom_MOE_1.fq in2=DR086_Erbom_MOE_2.fq out1=DR086_Erbom_MOE_1_corrected.fq out2=DR086_Erbom_MOE_2_corrected.fq -Xmx20g

Exception in thread "Thread-60" java.lang.NoClassDefFoundError: java/util/concurrent/ThreadLocalRandom

        at jgi.KmerNormalize$ProcessThread.run(KmerNormalize.java:2807)

Step 2: Deduplication

dedupe.sh in=NBS1170_Desrot_MOE_1_clean.fq out=NBS1170_Desrot_MOE_1_clean_dedupe.fq am ac fo renameclusters=f storename=f pto cc qin=33 csf=stats.txt dot=graph.dot -Xmx20g overwrite=t

dedupe.sh in=NBS1170_Desrot_MOE_2_clean.fq out=NBS1170_Desrot_MOE_2_clean_dedupe.fq am ac fo renameclusters=f storename=f pto cc qin=33 csf=NBS1170_Desrot_MOE_2_clean_dedupe_stats.txt dot=NBS1170_Desrot_MOE_2_clean_dedupe_graph.dot -Xmx20g overwrite=t

Step 3: Merge Paired End Reads

#merge PE 

bbmerge.sh in1=NBS1170_Desrot_MOE_1_clean.fq.gz in2=NBS1170_Desrot_MOE_2_clean.fq.gz out=NBS1170_Desrot_MOE_merged.fq -Xmx20g

Step 3.5: Could also remove the duplicates after merged paired ends

#remove duplicates of the two lanes

dedupe.sh in=NBS1170_Desrot_MOE_merged.fq out=NBS1170_Desrot_MOE_merged_dedupe.fq am ac renameclusters=f storename=f pto cc qin=33 csf=dedupe_logs/NBS1170_Desrot_MOE_merged_dedupe_stats.txt dot=dedupe_logs/NBS1170_Desrot_MOE_merged_dedupe_stats_graph.dot -Xmx20g overwrite=t

Compilation function of each BBmap shell scripts:

addadapters.sh

Randomly adds adapters to a file, or grades a trimmed file.

bbcountunique.sh

Generates a kmer uniqueness histogram, binned by file position.

There are 3 columns for single reads, 6 columns for paired:

count      number of reads or pairs processed

r1_first   percent unique 1st kmer of read 1

r1_rand    percent unique random kmer of read 1

r2_first   percent unique 1st kmer of read 2

r2_rand    percent unique random kmer of read 2

pair       percent unique concatenated kmer from read 1 and 2

bbduk.sh

Compares reads to the kmers in a reference dataset, optionally 

allowing an edit distance. Splits the reads into two outputs - those that 

match the reference, and those that don't. Can also trim (remove) the matching 

parts of the reads rather than binning the reads. Can use for quality trimming.

bbduk2.sh

Compares reads to the kmers in a reference dataset, optionally 

allowing an edit distance. Splits the reads into two outputs - those that 

match the reference, and those that don't. Can also trim (remove) the matching 

parts of the reads rather than binning the reads.

bbest.sh

Calculates EST (expressed sequence tags) capture by an assembly from a sam file.

Designed to use BBMap output generated with these flags: k=13 maxindel=100000 custom tag ordered

bbfakereads.sh

Generates fake read pairs from ends of contigs or single reads.

bbmap.sh

Fast and accurate short-read aligner for DNA and RNA.

bbmapskimmer.sh

Fast and accurate short-read aligner for DNA and RNA. (not sure difference from above)

bbmask.sh

Masks sequences of low-complexity, or containing repeat kmers, or covered by mapped reads.

bbmerge.sh

Merges paired reads into single reads by overlap detection.

With sufficient coverage, can also merge nonoverlapping reads using gapped kmers.

bbmergegapped.sh

Merges paired reads into single reads by overlap detection.

With sufficient coverage, can also merge nonoverlapping reads using gapped kmers.

bbnorm.sh

Normalizes read depth based on kmer counts.

Can also error-correct, bin reads by kmer depth, and generate a kmer depth histogram.

bbqc.sh

Performs quality-trimming; artifact, human, and phiX removal; adapter-trimming; error-correction and normalization. Designed for Illumina fragment libraries only.

bbrename.sh

Renames reads to <prefix>_<number> where you specify the prefix and the numbers are ordered.

bbsplit.sh

Maps reads to multiple references simultaneously.

Outputs reads to a file for the reference they best match, with multiple options for dealing with ambiguous mappings.

bbsplitpairs.sh

Separates paired reads into files of 'good' pairs and 'good' singletons by removing 'bad' reads that are shorter than a min length.

Designed to handle situations where reads become too short to be useful after trimming.  This program also optionally performs quality trimming.

bbwrap.sh

Wrapper for BBMap to allow multiple input and output files for the same reference.

calcmem.sh

????

calctruequality.sh

Calculates the observed quality scores from a sam file.

callpeaks.sh

No description--but appears related to keeping only certain reads within a certain region of histogram.

countbarcodes.sh

Counts the number of reads with each barcode.

countgc.sh

Counts GC content of reads or scaffolds.

crosscontaminate.sh

Generates synthetic cross-contaminated files from clean files.

Intended for use with synthetic reads generated by RandomReads.

cutprimers.sh

Cuts out sequences corresponding to primers identified in sam files.

decontaminate.sh

Decontaminates multiplexed assemblies via normalization and mapping.

dedupe.sh

Accepts one or more files containing sets of sequences (reads or scaffolds).

Removes duplicate sequences, which may be specified to be exact matches, subsequences, or sequences within some percent identity. Can also find overlapping sequences and group them into clusters.

dedupe2.sh

Accepts one or more files containing sets of sequences (reads or scaffolds).

Removes duplicate sequences, which may be specified to be exact matches, subsequences, or sequences within some percent identity.Can also find overlapping sequences and group them into clusters. (Not sure difference between dedupe.sh)

demuxbyname.sh

Demultiplexes reads based on their name (suffix or prefix) into multiple files.

ecc.sh

Corrects substitution errors in reads using kmer depth information.

Can also normalize and/or bin reads by kmer depth.

filterbarcodes.sh

Filters barcodes by quality, and generates quality histograms.

filterbycoverage.sh

Filters an assembly by contig coverage.

filterbyname.sh

Filters reads by name.

getreads.sh

Gets reads by number. The first read (or pair) has ID 0, the second read (or pair) has ID 1, etc.

grademerge.sh

Grades correctness of merging synthetic reads with headers generated by RandomReads and re-headered by RenameReads

gradesam.sh

Grades mapping correctness of a sam file of synthetic reads with headers generated by RandomReads3.java

idmatrix.sh

Generates an identity matrix via all-to-all alignment.

khist.sh

Generates a histogram of kmer counts for the input reads or assemblies.

Can also normalize, error-correct, and/or bin reads by kmer depth.

kmercount.sh

--counts kmers?

kmercountexact.sh

Counts the number of unique kmers in a file.

makechimeras.sh

Makes chimeric PacBio reads from nonchimeric reads.

mapPacBio.sh

Fast and accurate short-read aligner for DNA and RNA.

mapPacBio8k.sh

Description:  Fast and accurate short-read aligner for DNA and RNA.To index:   bbmap.sh ref=<reference fasta>

To map:     bbmap.sh in=<reads> out=<output sam>

To map without writing an index:

    bbmap.sh ref=<reference fasta> in=<reads> out=<output sam> nodisk

mapnt.sh

Maps sequences to the nt database.

matrixtocolumns.sh

Turns identity matrices into 2-column format for plotting.

mergeOTUs.sh

Merges coverage stats lines (from pileup) for the same OTU,

according to some custom naming scheme.

megebarcodes.sh

Concatenates barcodes and quality onto read names.

msa.sh

Aligns a query sequence to reference sequences.

Outputs the best matching position per reference sequence.

If there are multiple queries, only the best-matching query will be used.

phylip2fasta.sh

Calculates per-scaffold coverage information from an unsorted sam file.

pileup.sh

Calculates per-scaffold coverage information from an unsorted sam file.

printtime.sh

Prints time elapsed since last called on the same file.

randomreads.sh

Generates random synthetic reads from a reference genome.  Read names indicate their genomic origin. Allows precise customization of things like insert size and synthetic mutation type, sizes, and rates. Read names generated by this program are used by MakeRocCure (samtoroc.sh) and GradeSamFile (gradesam.sh). They can also be used by BBMap (bbmap.sh) and BBMerge (bbmerge.sh) to automatically calculate true and false positive rates, if the flag 'parsecustom' is used.

readlength.sh

Generates a length histogram of input reads.

reformat.sh

Reformats reads to change ASCII quality encoding, interleaving, file format, or compression format.

removehuman.sh

Removes all reads that map to the human genome with at least 95% identity after quality trimming.

removesmartbell.sh

Remove Smart Bell adapters from PacBio reads

repair.sh

Re-pairs reads that became disordered or had some mates eliminated. (not repair)

rqcfilter.sh

Performs quality-trimming, artifact removal, linker-trimming, adapter trimming, and spike-in removal using BBDukF. Performs human contaminant removal using BBMap.

samtoroc.sh

Creates a ROC curve from a sam file of synthetic reads with headers generated by RandomReads3.java

seal.sh

Performs high-speed alignment-free sequence quantification,

by counting the number of long kmers that match between a read and

a set of reference sequences.  Designed for RNA-seq with alternative splicing.

shuffle.sh

Reorders reads randomly.

stats.sh

Generates basic assembly statistics such as scaffold count, N50, L50, GC content, gap percent, etc.

staswraapper.sh

Runs stats.sh on multiple assemblies to produce one ouput line per file.

synthmda.sh

Generates synthetic reads following an MDA-amplified singe cell's coverage distribution.

testformat.sh

Tests the format of a sequence file based on name and contents.

textfile.sh

Translates nucleotide sequences to all 6 amino acid frames.

#installing Fasta_reader (needed to run script for stats on Trinity assembly)
#first need to install npm
brew install node
vi server.js
#paste the following
var http = require('http'); http.createServer(function (req, res) { res.writeHead(200, {'Content-Type': 'text/plain'}); res.end('Hello World\n'); }).listen(8124, "127.0.0.1"); console.log('Server running at http://127.0.0.1:8124/');
#ESC wq!

#to check if it works
node server.js

npm install fastareader

WOOPS! nevermind! All of this stuff is already installed with Trinity. To run TrinityStats.pl, it cannot be run alone. You must navigate to the directory in which it is found.

###################

#   Ar_jam MOE    #

###################

#combine data from Lanes 1 & 2

cat Arjam_MOE* >> DR_004_Arjam_MOE_assembled.fasta

sudo perl /Volumes/Spare/transcriptomes/Trinity/util/TrinityStats.pl /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/DR_004_Arjam_MOE_assembled.fasta

################################

## Counts of transcripts, etc.

################################

Total trinity transcripts: 140189

Total trinity components: 63956

Percent GC: 55.16

########################################

Stats based on ALL transcript contigs:

########################################

Contig N10: 4088

Contig N20: 3085

Contig N30: 2434

Contig N40: 1927

Contig N50: 1502

Median contig length: 426

Average contig: 826.50

Total assembled bases: 115865619

#####################################################

## Stats based on ONLY LONGEST ISOFORM per COMPONENT:

#####################################################

Contig N10: 3860

Contig N20: 2862

Contig N30: 2212

Contig N40: 1726

Contig N50: 1310

Median contig length: 440

Average contig: 793.55

Total assembled bases: 50752060 

###################

#   Ar_jam VNO    #

###################

#combine data from Lanes 1 & 2

cat Arjam_VNO* >> DR_011_Arjam_VNO_assembled.fasta

sudo perl /Volumes/Spare/transcriptomes/Trinity/util/TrinityStats.pl /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/DR_011_Arjam_VNO_assembled.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity transcripts: 138142
Total trinity components: 61477
Percent GC: 55.51

########################################
Stats based on ALL transcript contigs:
########################################

Contig N10: 4735
Contig N20: 3550
Contig N30: 2822
Contig N40: 2240
Contig N50: 1741

Median contig length: 455
Average contig: 920.11
Total assembled bases: 127106357

#####################################################
## Stats based on ONLY LONGEST ISOFORM per COMPONENT:
#####################################################

Contig N10: 3855
Contig N20: 2871
Contig N30: 2270
Contig N40: 1814
Contig N50: 1436

Median contig length: 467
Average contig: 848.36

Total assembled bases: 52154403

###################

#   Mo_red MOE    #

###################

#combine data from Lanes 1 & 2

cat Mored_MOE* >> DR_013_Mored_MOE_assembled.fasta

sudo perl /Volumes/Spare/transcriptomes/Trinity/util/TrinityStats.pl /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/DR_013_Mored_MOE_assembled.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity transcripts: 245573
Total trinity components: 100799
Percent GC: 56.14

########################################
Stats based on ALL transcript contigs:
########################################

Contig N10: 5311
Contig N20: 4036
Contig N30: 3304
Contig N40: 2703
Contig N50: 2118

Median contig length: 406
Average contig: 960.45
Total assembled bases: 235861464

#####################################################
## Stats based on ONLY LONGEST ISOFORM per COMPONENT:
#####################################################

Contig N10: 4140
Contig N20: 2880
Contig N30: 2105
Contig N40: 1512
Contig N50: 1053

Median contig length: 377
Average contig: 685.78

Total assembled bases: 69126334

###################

#   Mo_red VNO    #

###################

#combine data from Lanes 1 & 2

cat Mored_VNO* >> DR_013_Mored_VNO_assembled.fasta

sudo perl $TRINITY_HOME/util/TrinityStats.pl DR_013_Mored_VNO_assembled.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity transcripts: 163233
Total trinity components: 70319
Percent GC: 55.74

########################################
Stats based on ALL transcript contigs:
########################################

Contig N10: 4648
Contig N20: 3450
Contig N30: 2778
Contig N40: 2186
Contig N50: 1709

Median contig length: 430
Average contig: 883.37
Total assembled bases: 144195048

#####################################################
## Stats based on ONLY LONGEST ISOFORM per COMPONENT:
#####################################################

Contig N10: 3865
Contig N20: 2811
Contig N30: 2110
Contig N40: 1605
Contig N50: 1187

Median contig length: 427
Average contig: 754.88

Total assembled bases: 53082440

###################

#   Mo_bla MOE    #

###################

cat Mobl_MOE* >> DR_091_Mobl_MOE_assembled.fasta

sudo perl $TRINITY_HOME/util/TrinityStats.pl DR_091_Mobl_MOE_assembled.fasta
################################
## Counts of transcripts, etc.
################################
Total trinity transcripts: 249551
Total trinity components: 115871
Percent GC: 54.01

########################################
Stats based on ALL transcript contigs:
########################################

Contig N10: 4786
Contig N20: 3549
Contig N30: 2741
Contig N40: 2123
Contig N50: 1598

Median contig length: 335
Average contig: 763.67
Total assembled bases: 190573978

#####################################################
## Stats based on ONLY LONGEST ISOFORM per COMPONENT:
#####################################################

Contig N10: 3761
Contig N20: 2670
Contig N30: 1971
Contig N40: 1426
Contig N50: 989

Median contig length: 344
Average contig: 636.73

Total assembled bases: 73778325

###################

#Lane 2 Ar_jam MOE#

###################

gunzip pair*

cat pair.DR004_Artibeus_Jamaicensis_MOE_ACAGTG_L002_R1* >> Arjam_MOE_L2_R1.fq

cat pair.DR004_Artibeus_Jamaicensis_MOE_ACAGTG_L002_R2* >> Arjam_MOE_L2_R2.fq

fastq_quality_filter -i Arjam_MOE_L2_R1.fq -o Arjam_MOE_R1_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 15317036 reads.

Output: 14847480 reads.

discarded 469556 (3%) low-quality reads.

fastq_quality_filter -i Arjam_MOE_L2_R2.fq -o Arjam_MOE_R2_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 15317036 reads.

Output: 14608448 reads.

discarded 708588 (4%) low-quality reads.

~/Scripts/both.py Arjam_MOE_R1_L2_fastxtrimmed.fastq Arjam_MOE_R2_L2_fastxtrimmed.fastq

#run trinity

Trinity.pl --seqType fq --left Arjam_MOE_R1_L2_fastxtrimmed.fastq.both --right Arjam_MOE_R2_L2_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane2

###################

#Lane 1 Ar_jam VNO#

###################

gunzip pair*

cat pair.DR011_Artibeus_Jamaicensis_VNO_TGACCA_L001_R1* >> Arjam_VNO_L1_R1.fq

cat pair.DR011_Artibeus_Jamaicensis_VNO_TGACCA_L001_R2* >> Arjam_VNO_L1_R2.fq

mv Arjam_VNO_L1_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_011_Arjam_VNO/

fastq_quality_filter -i Arjam_VNO_L1_R1.fq -o Arjam_VNO_R1_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 13582991 reads.

Output: 13087556 reads.

discarded 495435 (3%) low-quality reads.

fastq_quality_filter -i Arjam_VNO_L1_R2.fq -o Arjam_VNO_R2_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 13582991 reads.

Output: 12665813 reads.

discarded 917178 (6%) low-quality reads.

~/Scripts/both.py Arjam_VNO_R1_L1_fastxtrimmed.fastq Arjam_VNO_R2_L1_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Arjam_VNO_R1_L1_fastxtrimmed.fastq.both --right Arjam_VNO_R2_L1_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane1

###################

#Lane 2 Ar_jam VNO#

###################

gunzip pair*

cat pair.DR011_Artibeus_Jamaicensis_VNO_TGACCA_L002_R1* >> Arjam_VNO_L2_R1.fq

cat pair.DR011_Artibeus_Jamaicensis_VNO_TGACCA_L002_R2* >> Arjam_VNO_L2_R2.fq

mv Arjam_VNO_L2_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_011_Arjam_VNO/

fastq_quality_filter -i Arjam_VNO_L2_R1.fq -o Arjam_VNO_R1_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 13746292 reads.

Output: 13242643 reads.

discarded 503649 (3%) low-quality reads.

fastq_quality_filter -i Arjam_VNO_L2_R2.fq -o Arjam_VNO_R2_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 13746292 reads.

Output: 12800383 reads.

discarded 945909 (6%) low-quality reads.

~/Scripts/both.py Arjam_VNO_R1_L2_fastxtrimmed.fastq Arjam_VNO_R2_L2_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Arjam_VNO_R1_L2_fastxtrimmed.fastq.both --right Arjam_VNO_R2_L2_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane2

###################

#Lane 1 Mo_red MOE#

###################

gunzip pair*

cat pair.DR013_Monophyllus_Redmani_MOE_CAGATC_L001_R1* >> Mored_MOE_L1_R1.fq

cat pair.DR013_Monophyllus_Redmani_MOE_CAGATC_L001_R2* >> Mored_MOE_L1_R2.fq

mv Mored_MOE_L1_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_013_Mored_MOE/

fastq_quality_filter -i Mored_MOE_L1_R1.fq -o Mored_MOE_R1_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 14546140 reads.

Output: 14205942 reads.

discarded 340198 (2%) low-quality reads.

fastq_quality_filter -i Mored_MOE_L1_R2.fq -o Mored_MOE_R2_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 14546140 reads.

Output: 13998979 reads.

discarded 547161 (3%) low-quality reads.

~/Scripts/both.py Mored_MOE_R1_L1_fastxtrimmed.fastq Mored_MOE_R2_L1_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mored_MOE_R1_L1_fastxtrimmed.fastq.both --right Mored_MOE_R2_L1_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane1

###################

#Lane 2 Mo_red MOE#

###################

gunzip pair*

cat pair.DR013_Monophyllus_Redmani_MOE_CAGATC_L002_R1* >> Mored_MOE_L2_R1.fq

cat pair.DR013_Monophyllus_Redmani_MOE_CAGATC_L002_R2* >> Mored_MOE_L2_R2.fq

mv Mored_MOE_L2_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_013_Mored_MOE/

fastq_quality_filter -i Mored_MOE_L2_R1.fq -o Mored_MOE_R1_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 14719223 reads.

Output: 14370971 reads.

discarded 348252 (2%) low-quality reads.

fastq_quality_filter -i Mored_MOE_L2_R2.fq -o Mored_MOE_R2_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v

Quality cut-off: 20

Minimum percentage: 80

Input: 14719223 reads.

Output: 14152391 reads.

discarded 566832 (3%) low-quality reads.

~/Scripts/both.py Mored_MOE_R1_L2_fastxtrimmed.fastq Mored_MOE_R2_L2_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mored_MOE_R1_L2_fastxtrimmed.fastq.both --right Mored_MOE_R2_L2_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane2

###################

#Lane 1 Mo_red VNO#

###################

gunzip pair*

cat pair.DR013_Monophyllus_Redmani_VNO_AGTCAA_L001_R1* >> Mored_VNO_L1_R1.fq

cat pair.DR013_Monophyllus_Redmani_VNO_AGTCAA_L001_R2* >> Mored_VNO_L1_R2.fq

mv Mored_VNO_L1_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_013_Mored_VNO/

fastq_quality_filter -i Mored_VNO_L1_R1.fq -o Mored_VNO_R1_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 14001874 reads.
Output: 13728472 reads.
discarded 273402 (1%) low-quality reads
fastq_quality_filter -i Mored_VNO_L1_R2.fq -o Mored_VNO_R2_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 14001874 reads.
Output: 13444029 reads.
discarded 557845 (3%) low-quality reads.

~/Scripts/both.py Mored_VNO_R1_L1_fastxtrimmed.fastq Mored_VNO_R2_L1_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mored_VNO_R1_L1_fastxtrimmed.fastq.both --right Mored_VNO_R2_L1_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane1

###################

#Lane 2 Mo_red VNO#

###################

gunzip pair*

cat pair.DR013_Monophyllus_Redmani_VNO_AGTCAA_L002_R1* >> Mored_VNO_L2_R1.fq

cat pair.DR013_Monophyllus_Redmani_VNO_AGTCAA_L002_R2* >> Mored_VNO_L2_R2.fq

mv Mored_VNO_L2_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_013_Mored_VNO/

fastq_quality_filter -i Mored_VNO_L2_R1.fq -o Mored_VNO_R1_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 14176309 reads.
Output: 13892340 reads.
discarded 283969 (2%) low-quality reads.
fastq_quality_filter -i Mored_VNO_L2_R2.fq -o Mored_VNO_R2_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 
-v
Quality cut-off: 20
Minimum percentage: 80
Input: 14176309 reads.
Output: 13594457 reads.
discarded 581852 (4%) low-quality reads.

~/Scripts/both.py Mored_VNO_R1_L2_fastxtrimmed.fastq Mored_VNO_R2_L2_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mored_VNO_R1_L2_fastxtrimmed.fastq.both --right Mored_VNO_R2_L2_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane2

###################

#Lane 1 Mo_bla MOE#

###################

gunzip pair*

cat pair.DR091_Mormoops_Blainvelli_MOE_CGATGT_L001_R1* >> Mobla_MOE_L1_R1.fq

cat pair.DR091_Mormoops_Blainvelli_MOE_CGATGT_L001_R2* >> Mobla_MOE_L1_R2.fq

mv Mobla_MOE_L1_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_091_Mobl_MOE/

fastq_quality_filter -i Mobla_MOE_L1_R1.fq -o Mobl_MOE_R1_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 16080335 reads.
Output: 15736163 reads.

discarded 344172 (2%) low-quality reads.
fastq_quality_filter -i Mobla_MOE_L1_R2.fq -o Mobl_MOE_R2_L1_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 16080335 reads.
Output: 15462533 reads.

discarded 617802 (3%) low-quality reads.

~/Scripts/both.py Mobl_MOE_R1_L1_fastxtrimmed.fastq Mobl_MOE_R2_L1_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mobl_MOE_R1_L1_fastxtrimmed.fastq.both --right Mobl_MOE_R2_L1_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane1

###################

#Lane 2 Mo_bla MOE#

###################

gunzip pair*

cat pair.DR091_Mormoops_Blainvelli_MOE_CGATGT_L002_R1* >> Mobla_MOE_L2_R1.fq

cat pair.DR091_Mormoops_Blainvelli_MOE_CGATGT_L002_R2* >> Mobla_MOE_L2_R2.fq

mv Mobla_MOE_L2_R* /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_091_Mobl_MOE/

fastq_quality_filter -i Mobla_MOE_L2_R1.fq -o Mobl_MOE_R1_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 16263150 reads.
Output: 15909002 reads.

discarded 354148 (2%) low-quality reads.
fastq_quality_filter -i Mobla_MOE_L2_R2.fq -o Mobl_MOE_R2_L2_fastxtrimmed.fastq -q 20 -p 80 -Q 33 -v
Quality cut-off: 20
Minimum percentage: 80
Input: 16263150 reads.
Output: 15621798 reads.

discarded 641352 (3%) low-quality reads.

~/Scripts/both.py Mobl_MOE_R1_L2_fastxtrimmed.fastq Mobl_MOE_R2_L2_fastxtrimmed.fastq

Trinity.pl --seqType fq --left Mobl_MOE_R1_L2_fastxtrimmed.fastq.both --right Mobl_MOE_R2_L2_fastxtrimmed.fastq.both --CPU 4 --JM 20G --output output_Lane2

Fail.

succeeded(0), failed(58736)   99.9949% completed.    libc++abi.dylib: terminating with uncaught exception of type std::length_error: vector::_M_fill_insert

succeeded(0), failed(58737)   99.9966% completed.    libc++abi.dylib: terminating with uncaught exception of type std::length_error: vector::_M_fill_insert

succeeded(0), failed(58738)   99.9983% completed.    libc++abi.dylib: terminating with uncaught exception of type std::length_error: vector::_M_fill_insert

succeeded(0), failed(58739)   100% completed.

We are sorry, commands in file: [failed_quantify_graph_commands.6513.txt] failed.  :-(

Error, cmd: /usr/local/Cellar/trinity/r20131110/trinity-plugins/parafly/bin/ParaFly -c /Volumes/Spare/transcriptomes/MOE_VNO_transcriptomes/Trinity_assembly/to_assemble/DR_004_Arjam_MOE/trinity_output//chrysalis/quantifyGraph_commands -CPU 4 -failed_cmds failed_quantify_graph_commands.6513.txt -v -shuffle  died with ret 256 at /usr/local/bin/Trinity.pl line 1793.

This is to be continued. At least they put a frowny face in the error message. 

Here is a blog that gives excellent Trinity instructions:
https://cartwrightlab.wikispaces.com/RNAseq+Assembly+in+Trinity

#for paired end reads, combine each side of pairs into one file

cat pair.DR004_Artibeus_Jamaicensis_MOE_ACAGTG_L001_R1* >> Arjam_MOE_R1.fq

cat pair.DR004_Artibeus_Jamaicensis_MOE_ACAGTG_L001_R2* >> Arjam_MOE_R2.fq

#run a trimming script

python ~/Scripts/q-trim.py Arjam_MOE_R1.fq Arjam_MOE_R1_trim.fq

#woops I run into this error:

importerror no module named screed

#it looks like I need the "screed" module for this; looks like I also need to update "setup tools"

#first install pip

sudo python get-pip.py

pip install ipython #not sure what this did

#still missing the screed module

sudo pip install setuptools --upgrade

sudo pip install git+https://github.com/ged-lab/screed.git

#okay seems to be working now

#trimmed with the default setting in which a minimum Qscore = 5

#minimum length = 30; lets see how this works

python ~/Scripts/q-trim.py Arjam_MOE_R1.fq Arjam_MOE_R1_trim.fq

python ~/Scripts/q-trim.py Arjam_MOE_R2.fq Arjam_MOE_R2_trim.fq

#now i need to assemble paired reads--extracting only the pairs using both.py

python ~/Scripts/both.py Arjam_MOE_R1_trim.fq Arjam_MOE_R2_trim.fq

#this creates two files ...R1_trim.fq.both and R2_trim.fq.both

#this script took about one hour to run

#finally time to run trinity!

Trinity.pl --seqType fq --left Arjam_MOE_R1_trim.fq.both --right Arjam_MOE_R2_trim.fq.both --CPU 4 --JM 20G --output trinity_output/

Everything is up and running now. Stay tuned.

Because I did such a horrid job taking comprehensible notes while in London, I am going to be more efficient this time around with transcriptome assembly.

We have our illumina RNA-seq reads back (yay!).

To download the reads from UArizona, you have to download this horridly complicated program called irods (http://uagc.arl.arizona.edu/resources/user-guide). After much struggle, the easiest way to get things installed and compiled is the "non-fancy" way. HOWEVER, you MUST download the most recent version of iRods (as of right now 3.3.1)

Navigate to your iRODs folder. Things in blue are what you type. Things in pink are what you would change differently based on your computer.

 export PATH=$PATH:$HOME/iRODS/clients/icommands/bin

 iinit #this will prompt you for a bunch of information

Enter the host name (DNS) of the server to connect to:  mozart.arl.arizona.edu

Enter the port number:  1247

Enter your irods user name:  <your login name from UAGC/BCF>

Enter your irods zone:  BCFZone

Enter your current iRODS password:  <your password from UAGC/BCF>

 

#now your session should start

ils #lists all of the folders in your iRODS account

icd 14140709_SN885_0197_Lane1SideA_Project_Davalos_RNA #change directory

iget -r 140709_SN885_0197_Lane1SideA_Project_Davalos_RNA /Volumes/Spare/ 

#iget -r copies the folder you want to whatever directory you want

Okay now all of the data should be on your computer.

Now I need to install Trinity:

#install homebrew if you do not already have it

ruby -e "$(curl -fsSL https://raw.github.com/Homebrew/homebrew/go/install)"

#now install programs necessary for trinity to run on your mac

brew install homebrew/science/bowtie

brew install git

#to allow for new installation of compiler

brew tap homebrew/dupes

brew install apple-gcc42

#okay guess I can't do that because it requires an operating system higher than Mac OSX Lion...

Yes, you read that correctly. Our operating system is still pre-lion. Seems like I will have to wait until our operating system is figured out.

Once we figure this out though, remember this link: http://sourceforge.net/p/trinityrnaseq/mailman/message/31892980/

It allows for trinity to be installed on Mac (even though it is compiled for Linux).