ECCB 2016 Den Hague, Netherlands computational Biologist and Bioinformaticians gatherings at a sweet Dutch country !

I have been to several conferences within India, while ECCB 2016 which happened in Den Hague, from September 3 2016 - September 7 2016. It was the first time for me travel outside India, had butterflies on my stomach the day before I travel.The trip really went well. It was a gathering of computational Biologist and Bioinformaticians over the world. Well I should thank Department of Science and Technology, Government of India for providing me the travel award. The Meeting started with the workshop on discussing Pacbio and Nanopore data. Expertise from the field of nanopore and Pacbio were discussing the problems with the long reads. People were complaining about the "error rates" of these reads, and difficulties in genome assembly of these reads. Had a great opportunity to discuss with the experts. The Nanopore experts were suggesting that Canu assembler can do better when handling the problematic regions in the genome. The Miniasm and Racon assembler also be tried . There were sessions about Irys to create a genome map and align the created  map back to the genome assembly to get a better genome assembly. The structural variants and the comparative genomics are also studied from the graph. Next topic was using Isoseq from pacific Bio-systems to produce a full length transcripts without assembly, followed by promethion and squiggle sequencing system from nanopore technology"Read Until " approach it enables selection of individual DNA molecules for sequencing from a pool of DNA molecules. Then there was a session of Minotour where the base calling of the nanopore reads where done without performing the cloud base calling since there is a dependency of high speed internet. The developers of the tools and technologies were very friendly and gave suggestions on working with the long reads. after the workshop the conference scientific sessions started and many interesting talks where there, I was more interested towards the error correction algorithm development, genome assembly tools, new ortholog prediction tools. Most of the sessions and posters where about the cancer ( a devil), and ENCODE. I can say that 60% of people presented towards cancer transcriptomics and genomics, and 30 % of work in ENCODE, rest where like plant, bacteria, database development, Docking and simulations. The talks and discussions can be retrieved from twitter via #ECCB2016.  I liked the theme of the conference here not only PhD students and Scientists were presenting the work, even people working in companies were also showcased the ongoing work.Some people were very happy and showed interest towards the  poster of my PhD work, since its a plant pathogen. I am more interested towards studying the environmental organisms, pathogens of human, discovering various new species from the environment. About the food it was good, had a varieties of cheese. I had time to visit Amsterdam its a very nice place with a polite people. Visited churches, Museum, had a good canal Boat riding.I had few friends from conference and joined with them and rented a boat and rode over the Amsterdam city. The future ECCB2017 conference will be held in Prague.

Analyzing Differential expression analysis data using the tuxedo suite (cummeRbund)

Tuxedo suite comprises of bowtie, tophat, cufflink, cummeRBund and many more accessory tools.

First get your genome fasta file (final genome assembly file).
1. Map your RNAseq fastq files using tophat (if all is well your run will be seamless)
2. Run cufflink over your tophat output file (cufflinks accepted_hits.bam). This run will take a while since cufflink will actually merge the reads into transcripts, isoforms, genes and so on. If your files are large then in a good enough server expect it to run for 8-12 hours.
3. Run cuffmerge: cuffmerge list.txt -> where list.txt carries the names of the files of *_transcripts.gtf files. This will run very fast and will merge all the gene_ids that will be same across all your samples. The output of this file is a merged.gtf file.
4. For running differential expression analysis run the following:
/cuffdiff merged.gtf tophat_HTI1-vs-HTI4/accepted_hits.bam tophat_HTI2-vs-HTI4/accepted_hits.bam tophat_HTI3-vs-HTI4/accepted_hits.bam

This will create a plethora of files, but the following files are the ones you will be proceeding with for cummeRbund for result visualization and generating publication quality images.

For running cummeRbund, get all these files to your working directory

isoforms.fpkm_tracking

isoform_exp.diff

genes.fpkm_tracking

gene_exp.diff

tss_groups.fpkm_tracking

tss_group_exp.diff

cds.fpkm_tracking

cds_exp.diff

cds.diff

promoters.diff

splicing.diff

The best option will be to put all of these 11 files into a separate directory inside your working directory: say 'diff_exp'
You can run Rstudio if you like in your windows machine or run R in your server. For running CummeRbund you will need the following packages that you can go ahead and download upfront:

• RSQLite
• ggplot2 v0.9.2
• reshape2
• plyr
• fastcluster
• rtracklayer
• Gviz
• BiocGenerics (>=0.3.2)
• Hmisc

In case you have forgotten how to install R packages go this way: source('http://www.bioconductor.org/biocLite.R') biocLite('cummeRbund') And follow this same protocol for installing other R packages. Once done you can start with setting your working directory using setwd() command.
 For example: setwd("C:/Users/Sucheta/Documents/MyLabIICB/AllCollaborations/NahidAliCollaboration/companion")
Then load the library:

library(cummeRbund)

Now read your 11 files using this command
data <-readCufflinks("diff_exp")

This will take a while to read but will create a db file in your source directory. This is your database file.

Now you can plot gene density using the following command:

csDensity(genes(data))

Or can do a volcano plot of differentially expressed genes using:

v<-csVolcanoMatrix(genes(data))
v

As you can see from this file, the different conditions have least difference among themselves.

This will continue in next blog...

#OMGN2016 Malmo, Sweden - Between Then and Now...

Many things have changed in the years in front of me since the day I started attending OMGN meetings. My first meeting was in year 2005 and then the first Oomycetes genomes were getting sequenced and getting analyzed - at its own pace (read very slow pace). We used to get excited even when we got SSRs or repeats predicted. I distinctly remember the 2004 Joint Genome Institute sequence jamboree when in the evenings we used to gather to discuss what was done during the day. On second or third day of Jamboree, Brett came up with this multiple sequence alignment that presumably indicated that there was an RXLR motif in the effector proteins. It was a huge deal then. Subsequently in all the meetings everybody started discussing on these proteins. Initially it appeared too good to be true with this small 4 letter motif, but a lot of work was done especially in Brett's lab to prove that it indeed was a significant motif. The prediction algorithms of effectors got published in high flying journals, everybody was excited. Slowly many more papers came out on RXLRs, their prediction methods, characterizations till 2010. Now in 2016, I see the level of science has gone up way higher. Genome sequencing using PacBio or Illumina is no deal, neither is analyzing them. Effector prediction has become just a days job (Thanks to all the hard work done by the pioneers). Genome analysis are now carried out by single individuals in few months time. This meeting was a skew towards miRNAs, CRISPER technology for genome editing, pacbio sequencing, RNA silencing. The effector biology has moved many many steps up now. Many more things are now known. Many more proteins have been characterized. It is exciting time in the history of oomycetes biology where many things are happening right in front of our eyes. For those who could not attend please check #OMGN16 for more details. For me now bye bye lovely Malmo!!

Variant Calling - The bowtie - picard - samtools - gatk pipeline....

Nextgen sequencing has caused a sudden surge in data deluge, but the informatics pipelines and algorithms are unable to keep up with the pace. While most of the exome sequencing data finally focuses on SNP calling and there are various ways of doing this, I decided to discuss one pipeline that has been accepted all over as one of the most sophisticated methods. It is the bowtie - picard - gatk pipeline.

When you are dealing with colorspace data the choice of mappers get limited. Howver, my favorite mapper is still bowtie for several reasons. Lifescope has its own inhouse mapper; which claims to have a all round better approach in mapping colorspace data, but the lack of transparency on what happens within puts me off using this tool. Once bowtie maps the reads by default parameter, the next thing to do is to convert the sam file into bam file, sort it and index it. All these can be done using samtools. However, if the file size is large, you could do the sorting job using sort operations from unix commandline.

1. sort sam file
export TMPDIR=DIR_WITH_LOTS_OF_SPACE
LC_ALL="C" sort -k 3,3 -k 4,4n input_sam > output_sam # This step will take a long time

sort options for samtools works but only on bam files and on many instances downstream analysis softwares complain about co-ordinates not being sorted...

or Use picard:

java -jar /share/apps/picard-tools-1.56/SortSam.jar I=bowtie.sam O=bowtie.bam SO=coordinate # This took one hour in a HPC with 48 GB RAM on each node for a file size of 30 GB

2. Make an index file of bam file
samtools index bowtie.bam bowtie.bai

3. MarkDuplicates using picard
java -jar /share/apps/picard-tools-1.56/MarkDuplicates.jar I=bowtie.bam M=metrics.bam O=duplicateMarked.bam

4. sort this bam file and make index using samtools
samtools sort duplicateMarked.bam duplicateMarked.sorted
samtools index duplicateMarked.sorted.bam duplicateMarked.sorted.bai

5. Then run IndelRealignerTargetCreator using GATK
java -jar /share/apps/GenomeAnalysisTK-2.4-9-g532efad/resources/GenomeAnalysisTK.jar -T RealignerTargetCreator -I duplicateMarked.sorted.bam  -R /share/reference/human/samtools/hg19.fa -o dM.bam.list 
# The output file dM.bam.list returns 0 output. Check it later.

6. Now run this picard tool to get RG updated since indelRealigner complains.

java -jar /share/apps/picard-tools-1.56/AddOrReplaceReadGroups.jar I=duplicateMarked.bam O=readGroupReplaced.bam RGLB="LINK_TO_FASTA" RGPL=SOLID RGPU=run barcode RGSM=9111 SORT_ORDER=coordinate CREATE_INDEX=TRUE VALIDATION_STRINGENCY=LENIENT'

Posters I could take pictures of Beyond Genome 2014

Here are few of the posters in Beyond Genome meeting that I could take pictures of. There were many more, but access to take their pictures was less...

Our Poster

  

SNP calling using GATK for de novo genome

I have a got a chance to work in leishmania genome, where i have a genome assembly and i dont have any deposited dbSNP or any other reference file to do variant calling, i have been working and stuck in many steps and posted in GATK forums they replied to some of my queries at one point stopped to reply since People were having a fat and busy holiday on thanks giving , and figured out how to do the variant calling, i think this blog will be much useful for the naive person like me, lets see the workflow and please refer GATK documentation for the explanation.
#first build the index for the reference genome
/share/apps/bowtie2-2.1.0/bowtie2-build after_removing_2k.fasta  leishmania.index.bt2
#after index map the reference to the reads
/share/apps/bowtie2-2.1.0/bowtie2 -x leishmania.index.bt2 -1 /data/results/STLab/NahidAli/141218_SND393_A_L005_HTI-5_trim_R1_filtered.fastq -2 /data/results/STLab/NahidAli/141218_SND393_A_L005_HTI-5_trim_R2_filtered.fastq -S bowtie_aligned.sam
#convert the sam file to bam
samtools view -S bowtie_aligned.sam -b -o bowtie_aligned.bam
#sort the bam file
samtools sort bowtie_aligned.bam bowtie_aligned_sorted
#create a pileup file
samtools mpileup -uf after_removing_2k.fasta bowtie_aligned_sorted.bam|/share/apps/samtools-0.1.18/bcftools/bcftools view -bvcg - > leishmania.raw.bcf
#convert bcf to vcf
/share/apps/samtools-0.1.18/bcftools/bcftools view leishmania.raw.bcf > leishmania.raw.vcf
*********************************************************************************

 The bove commands are just initial way of mapping the reads to the reference and the real GATK pipeline starts below since i don't have the any known sites  i have done without base cailbration 

#create the dictionary
#java -jar /share/apps/picard-tools-1.56/CreateSequenceDictionary.jar R=after_removing_2k.fasta O=after_removing_2k.dict
#add or mark group ids
#java -jar /share/apps/picard-tools-1.56/AddOrReplaceReadGroups.jar I=bowtie_aligned.sam O=group_added_read.bam SO=coordinate RGID=1 RGLB=library1 RGPL=illumina RGPU=1 RGSM=leishmania VALIDATION_STRINGENCY=LENIENT CREATE_INDEX=TRUE
#mark duplicates
#java -jar /share/apps/picard-tools-1.56/MarkDuplicates.jar I=group_added_read.bam O=mapped_reads_dup.bam METRICS_FILE=metricsFile CREATE_INDEX=true
#sort bam file
#java -jar /share/apps/picard-tools-1.56/BuildBamIndex.jar INPUT=mapped_reads_dup.bam
#create realign target creator
#share/apps/GenomeAnalysisTK.jar -T RealignerTargetCreator -R after_removing_2k.fasta -o target_interval.intervals -I mapped_reads_dup.bam
#indel realigner
#/share/apps/GenomeAnalysisTK.jar -T IndelRealigner -R after_removing_2k.fasta -I mapped_reads_dup.bam -targetIntervals target_interval.intervals -o Indel_realigned.bam
#haplotype caller
#java -jar /share/apps/GenomeAnalysisTK.jar -T HaplotypeCaller -R after_removing_2k.fasta -I Indel_realigned.bam -stand_call_conf 30 -stand_emit_conf 10 -o raw_variants.vcf
#choose the variants from the raw vcf file
#java -jar /share/apps/GenomeAnalysisTK.jar -T SelectVariants -R after_removing_2k.fasta -V raw_variants.vcf -selectType SNP -o raw_snps.vcf
#do the filtration
#java -jar /share/apps/GenomeAnalysisTK.jar -T VariantFiltration -R after_removing_2k.fasta -V raw_variants.vcf --filterExpression "QD < 2.0 || FS > 60.0 || MQ < 40.0 || MQRankSum < -12.5 || ReadPosRankSum < -8.0" --filterName "my_snp_filter" -o filtered_snps.vcf
#Extract the indels
#java -jar /share/apps/GenomeAnalysisTK.jar -T SelectVariants -R after_removing_2k.fasta -V raw_variants.vcf -selectType INDEL -o raw_indels.vcf
#do the filteration
#java -jar /share/apps/GenomeAnalysisTK.jar -T VariantFiltration -R after_removing_2k.fasta -V raw_variants.vcf --filterExpression "QD < 2.0 || FS > 200.0 || ReadPosRankSum < -20.0" --filterName "my_indel_filter" -o filtered_indels.vcf

from the above predicted snps and indels extract the regions and further annotate and work on it happy variant calling !!!!!!!!!!

5C bed file data format

5C and 3C are the newer technologies in sequencing where the chromatin inetraction data can be obtained. If you looking for such data and happen to download from UCSC genome browser, it may be hard to look around for format describing the fields. We asked the authors and here is the explanation:

The site from which you may download data may be this: https://www.encodeproject.org/experiments/ENCSR000CYD/

BED  file format descrition can be found from : https://genome.ucsc.edu/FAQ/FAQformat.html#format1 

Here is a sample data for GM12878 cell line:

chr22   31998728        33247041        5C_301_ENm004_FOR_292.5C_301_ENm004_REV_

32      1000    .       31998728        33247041        0       2       12744,40

98,     0,1244215,

chr5    131346229       132145236       5C_299_ENm002_FOR_241.5C_299_ENm002_REV_

33      1000    .       131346229       132145236       0       2       2609,210

5,      0,796902,

col1: Chromosome name

col2: Chromosome start

col3: chromosome end

col4: Name of the interacting sites (primer names)

col5:

col7: chromosome start

col8: chromosome end

col11: block sizes in comma separated list

col12: block offset in comma separated list

Now I will explain what col11 and col12 means...

the beginning of interacting site is the cromosome start and the beginning of offset is 0.

So, the interacting site begins at 31998728 + 0 and the interacting block length is 12744.

The beginning position of interacting site 2 is: 31998728 + 1244215 = 33242943

 The size of interacting block 2 is 4098. so, end of interacting site is 33242943 + 4098 = 33247041.

Here is a diagrammatic representation: