18 September 2012

Notes about bwa 0.6.2 and the multiple hits.

Here are my notes about the way bwa 0.6.2 handles the multiple hits.


I've created a sub-sequence (~6Mb ) of the human chr22 (see the Makefile below). The sequence doesn't contain any base 'N' or any lowercase (=repeatMasked) letter.
1E6 reads (forward and reverse) have been generated from this ~chr22 with samtools wgsim. No mutation and no sequencing error have been allowed.

This ~chr22 sequence have been duplicated and renamed as 'DUP'. Both sequences have been merged in one fasta file named 'twoChrom.fa'

The reads have been aligned with BWA (v. 0.6.2-r126) on this genome.


There's only one SAM alignment per read

$ samtools view align.sorted.bam  | wc -l
2000000

1% of the reads were not properly paired/mapped on the genome

$ samtools view  -f 2 align.sorted.bam | wc -l
1999926

The properly mapped reads contain some informations about the alternative hits in 'XA'

$samtools view  -f 2 align.sorted.bam | grep 22_10_402_0 | verticalize -n

>>> 1
$1   22_10_402_0:0:0_0:0:0_c1608
$2   99
$3   22
$4   10
$5   0
$6   70M
$7   =
$8   333
$9   393
$10  GGGACGGTCATGCAATCTGGACAACATTCACCTTTAAAAGTTTATTGATCTTTTGTGACATGCACGTGGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:DUP,+10,70M,0;
<<< 1

>>> 2
$1   22_10_402_0:0:0_0:0:0_c1608
$2   147
$3   22
$4   333
$5   0
$6   70M
$7   =
$8   10
$9   -393
$10  ATACCTGCCAGATGAGTCACTGGCAAAAGGTGCTGCTCCCTGGTGAGGGAGAAACACCAGGGGCTGGGAG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:DUP,-333,70M,0;
<<< 2

50% mapped on chr22, 50% mapped on DUP

$ samtools view  -f 2 align.sorted.bam | cut -d '      ' -f 3 | sort | uniq -c
 999860 22
1000066 DUP

Some pairs have been unmapped because the mate have been mapped on the other chromosome !

$ samtools-0.1.18/samtools view  -F 2 align.sorted.bam | grep 22_1210321_1210924 | verticalize -n

>>> 1
$1   22_1210321_1210924_0:0:0_0:0:0_e28f5
$2   81
$3   22
$4   1210855
$5   0
$6   70M
$7   DUP
$8   1210321
$9   0
$10  AGAGACTGAGTCCGGCTAGAGAACAGGGTGGAGCCCCTTTGGACCTTAGAGCTGGGCCTTTGGGCCTTGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:6
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:22,-1996342,70M,0;DUP,-1996342,70M,0;22,+2551778,70M,0;DUP,+2551778,70M,0;DUP,-1210855,70M,0;
<<< 1

>>> 2
$1   22_1210321_1210924_0:0:0_0:0:0_e28f5
$2   161
$3   DUP
$4   1210321
$5   0
$6   70M
$7   22
$8   1210855
$9   0
$10  CCCGGGCCGGATGGCTCGCCTGCGCGGCCAGCTCCGGGCCGAAGCGGCTTCGCGGTCCGAGGTGCCGCGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:22,+1210321,70M,0;
<<< 2

The Makefile

SAMDIR=/usr/local/package/samtools-0.1.18
SAMTOOLS=$(SAMDIR)/samtools
BCFTOOLS=$(SAMDIR)/bcftools/bcftools
BWA=/usr/local/package/bwa-0.6.2/bwa
REF1=chr22.fa
REF2=twoChrom.fa

.INTERMEDIATE : align.sam random_1.sai random_2.sai align.bam variations.bcf

%.bam : %.sam
        $(SAMTOOLS) view -o $@ -b -S -T $(REF2) $<
%.bam.bai : %.bam
        $(SAMTOOLS) index $<



align.sorted.bam : align.bam
        $(SAMTOOLS) sort $< align.sorted


align.sam : random_1.sai random_2.sai  
        $(BWA) sampe -a 600 $(REF2) $^ random_1.fq.gz random_2.fq.gz > $@

$(REF1):
        curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/$(REF1).gz" |\
        gunzip -c | sed 's/^>chr/>/' | grep -v '[Nnatgc]' | head -n 100000 > $@


random_1.sai :  random_1.fq.gz $(REF2).bwt
        $(BWA) aln -f $@ $(REF2) $<

random_2.sai :  random_2.fq.gz $(REF2).bwt
        $(BWA) aln -f $@ $(REF2) $<

random_1.fq.gz random_2.fq.gz : $(REF1)
        $(SAMDIR)/misc/wgsim -N 1000000 -e 0.0 -r 0.0 -d 400 $< random_1.fq random_2.fq > wgsim.output
        gzip -f --best random_1.fq random_2.fq

$(REF2): $(REF1)
        cp $< $@
        sed 's/^>.*/>DUP/' $< >> $@

$(REF2).bwt : $(REF2)
        $(BWA) index -a bwtsw $<

$(REF2).fai :  $(REF2)
        $(SAMTOOLS) faidx $< 



clean:
        rm -f chr22.* *.bam *.vcf *.bcf *.sai *.gz *.fq *.bai  wgsim.output *.sam
That's it,

Pierre

13 September 2012

Translating a DNA sequence in Google Spreadsheet using #GoogleAppScript

Google has released Google App Script.
"Google Apps Script is a JavaScript cloud scripting language that lets you extend Google Apps and build web applications. Scripts are developed in Google Apps Script’s browser-based script editor, and they are stored in and run from Google's servers.
Google Apps Script is very versatile. Here are some examples of things you can do with Google Apps Script":
  • Build custom functions in a Google Spreadsheet
  • Extend certain Google Apps products by creating custom menus linked to scripts
  • Create and publish web applications, which can run on their own or embedded within a Google Site
  • Schedule tasks like report creation and distribution and run them on a custom schedule
  • Automate workflows such as document or expense approvals, order fulfillment, time-tracking, and more

In the current post, I'll show how to create a custom javascript function that will

Translate a
DNA
to a
Protein
into a
Google Spreadsheet

Create a new Google Spreadhseet. Open the menu "Tools" > "Script Manager...". Click on New...

Click on "Create Blank Project".

A new editor is opened. Copy the following javascript code (https://gist.github.com/3716137) into the editor. Save the javascript projet.

Close the script, go back to the spreasheet. You can now use your new function =translateDNA(dna):

That's it,

Pierre

05 September 2012

Customizing the java classes for the NCBI generated by XJC

Reminder: XJC is the Java XML Binding Compiler. It automates the mapping between XML documents and Java objects:

XSD (aka XML-schema)
+
XJC
=
JAVA Classes

The code generated by XJC allows to :
  • Unmarshal XML content into a Java representation
  • Access and update the Java representation
  • Marshal the Java representation of the XML content into XML content

For example, the following XML-Schema (tinyseq.xsd) describes a TinySeq-XML document returned by the NCBI.

<?xml version="1.0" encoding="UTF-8"?>
<xs:schema xmlns:xs="http://www.w3.org/2001/XMLSchema">
  <xs:annotation>
    <xs:documentation> XML schema for NCBI tinyseq format</xs:documentation>
  </xs:annotation>
  <xs:complexType name="TSeqSet_t">
    <xs:annotation>
      <xs:documentation>Set of sequences</xs:documentation>
    </xs:annotation>
    <xs:sequence>
      <xs:element ref="TSeq" maxOccurs="unbounded"/>
    </xs:sequence>
  </xs:complexType>
  <xs:complexType name="TSeq_t">
    <xs:annotation>
      <xs:documentation>A Tiny Sequence</xs:documentation>
    </xs:annotation>
    <xs:sequence>
      <xs:element name="TSeq_seqtype">
        <xs:complexType>
          <xs:attribute name="value">
            <xs:simpleType>
              <xs:restriction base="xs:string">
                <xs:enumeration value="nucleotide"/>
                <xs:enumeration value="protein"/>
              </xs:restriction>
            </xs:simpleType>
          </xs:attribute>
        </xs:complexType>
      </xs:element>
      <xs:element name="TSeq_gi" type="xs:long"/>
      <xs:element name="TSeq_accver" type="xs:string"/>
      <xs:element name="TSeq_sid" type="xs:string"/>
      <xs:element name="TSeq_taxid" type="xs:long"/>
      <xs:element name="TSeq_orgname" type="xs:string"/>
      <xs:element name="TSeq_defline" type="xs:string"/>
      <xs:element name="TSeq_length" type="xs:nonNegativeInteger"/>
      <xs:element name="TSeq_sequence" type="xs:string"/>
    </xs:sequence>
  </xs:complexType>
  <xs:element name="TSeqSet" type="TSeqSet_t"/>
  <xs:element name="TSeq" type="TSeq_t"/>
</xs:schema>

This xml-schema can be compiled with XJC:

${JAVA_HOME}/bin/xjc -d . -p generated tinyseq.xsd
parsing a schema...
compiling a schema...
generated/ObjectFactory.java
generated/TSeqSetT.java
generated/TSeqT.java

$ more generated/TSeqT.java

package generated;
(...)
@XmlAccessorType(XmlAccessType.FIELD)
@XmlType(name = "TSeq_t", propOrder = {  "tSeqSeqtype",  "tSeqGi",  "tSeqAccver","tSeqSid",(...),"tSeqSequence"})
public class TSeqT {
    @XmlElement(name = "TSeq_seqtype", required = true)
    protected TSeqT.TSeqSeqtype tSeqSeqtype;
    @XmlElement(name = "TSeq_gi")
    protected long tSeqGi;
    @XmlElement(name = "TSeq_accver", required = true)
    protected String tSeqAccver;
    @XmlElement(name = "TSeq_sid", required = true)
    protected String tSeqSid;
    @XmlElement(name = "TSeq_taxid")
    protected long tSeqTaxid;
    @XmlElement(name = "TSeq_orgname", required = true)
    protected String tSeqOrgname;
    @XmlElement(name = "TSeq_defline", required = true)
    protected String tSeqDefline;
    @XmlElement(name = "TSeq_length", required = true)
    @XmlSchemaType(name = "nonNegativeInteger")
    protected BigInteger tSeqLength;
    @XmlElement(name = "TSeq_sequence", required = true)
    protected String tSeqSequence;
(...)
    }

But XJC doesn't know how to generate some classical java functions like 'hashCode', 'equals' or 'toString' or to add some custom methods to your classes.

Hopefully the standard distribution of XJC comes with a plugin named -Xinject-code whch injects some custom code in the classes generated by XJC.

XSD (aka XML-schema)
+
java xml binding (jxb)
+
XJC
=
Customized JAVA Classes

For example, if we want to add a toString method to the class TSeqT, we're going to write the following "java xml binding file" (jxb) which alters the initial xml schema:

<?xml version="1.0" encoding="UTF-8"?>
<jxb:bindings 
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xmlns:xjc="http://java.sun.com/xml/ns/jaxb/xjc"
xmlns:jxb="http://java.sun.com/xml/ns/jaxb"
xmlns:ci="http://jaxb.dev.java.net/plugin/code-injector"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
jxb:extensionBindingPrefixes="ci "
jxb:version="2.1"

>

<jxb:bindings schemaLocation="tinyseq.xsd">
        <!-- here we use an XPATH expression to tell xjc about which part
 of the XML schema we want to change -->
 <jxb:bindings node="/xs:schema/xs:complexType[@name='TSeq_t']">
  <ci:code>

 /** toString : returns the gi and the defline  */
 public String toString()
  {
  return "gi:"+getTSeqGi()+"|"+getTSeqDefline();
  }

</ci:code>
 </jxb:bindings>
</jxb:bindings>

</jxb:bindings>

Below, I wrote a larger JXB file 'tinyseq.jxb' which injects the following methods:
  • 'equals' method for TSeq
  • 'hashCode' method for TSeq
  • 'toString' method for TSeq
  • 'printAsFasta' method for TSeq
  • 'getTSeqSetbyId' method for TSeqSet. A static function fetching a TinySeq sequence from the NCBI for a given 'gi'
  • a 'main' method for TSeqSet. It loops over a list of 'gi's, fetches the sequences (using NCBI-EFetch) and prints the sequences as FASTA
<?xml version="1.0" encoding="UTF-8"?>
<jxb:bindings 
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xmlns:xjc="http://java.sun.com/xml/ns/jaxb/xjc"
xmlns:jxb="http://java.sun.com/xml/ns/jaxb"
xmlns:ci="http://jaxb.dev.java.net/plugin/code-injector"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
jxb:extensionBindingPrefixes="ci "
jxb:version="2.1"

>

<jxb:bindings schemaLocation="tinyseq.xsd">
 <jxb:bindings node="/xs:schema/xs:complexType[@name='TSeq_t']">
  <ci:code>
 /* print this sequence as fasta */
 public void printAsFasta(java.io.PrintStream out)
  {
  String s=getTSeqSequence();
  out.print(">gi:"+getTSeqGi()+"|"+ getTSeqAccver() +"|"+getTSeqDefline());
  for(int i=0;i &lt; s.length();++i)
   {
   if(i%60==0) out.println();
   out.print(s.charAt(i));
   }
  out.println();
  }

 /** equals: two  TSeq are equal if they have the same gi */
 @Override
 public boolean equals(Object o)
  {
  if(o==this) return true;
  if(o==null || o.getClass()!=this.getClass()) return false;
  return this.getTSeqGi()==TSeqT.class.cast(o).getTSeqGi();
  }

 /** hashCode : use gi */
 @Override
 public int hashCode()
  {
  return  (int)(this.getTSeqGi()^(this.getTSeqGi()>>>32));
  }

 /** toString : returns the gi and the defline  */
 public String toString()
  {
  return "gi:"+getTSeqGi()+"|"+getTSeqDefline();
  }

</ci:code>
 </jxb:bindings>




 <jxb:bindings node="/xs:schema/xs:complexType[@name='TSeqSet_t']">
  <ci:code>
 /** get TSeqSetT from a given gi */
 public static TSeqSetT getTSeqSetbyId(long gi)
  throws javax.xml.bind.JAXBException, javax.xml.bind.UnmarshalException , java.io.IOException
  {
  /** find the JAXB context in the defined path */
  javax.xml.bind.JAXBContext jc = javax.xml.bind.JAXBContext.newInstance(TSeqSetT.class,TSeqT.class);
  javax.xml.bind.Unmarshaller u = jc.createUnmarshaller();
  String uri="http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&amp;rettype=fasta&amp;retmode=xml&amp;id="+gi;
  /** read the sequence */
  return u.unmarshal(new javax.xml.transform.stream.StreamSource(uri),TSeqSetT.class).getValue();
  }

 /** main: takes a list of gi and prints the sequences as fasta */
 public static void main(String args[]) throws Exception
  {
  for(int optind=0;optind &lt; args.length;++optind)
   {
   TSeqSetT tss=TSeqSetT.getTSeqSetbyId(Long.parseLong(args[optind]));
   for(generated.TSeqT seq:tss.getTSeq()) seq.printAsFasta(System.out);
   }
  }

</ci:code>

</jxb:bindings>

</jxb:bindings>

</jxb:bindings>

Compile the schema, compile the java classes and execute

$ xjc -target 2.1 -verbose -Xinject-code -extension -d . -p generated -b tinyseq.jxb tinyseq.xsd
parsing a schema...
compiling a schema...
[INFO] generating code
unknown location

generated/ObjectFactory.java
generated/TSeqSetT.java
generated/TSeqT.java

$ javac generated/*.java
$ java generated.TSeqSetT 25 26 27
>gi:25|X53813.1|Blue Whale heavy satellite DNA
TAGTTATTCAACCTATCCCACTCTCTAGATACCCCTTAGCACGTAAAGGAATATTATTTG
GGGGTCCAGCCATGGAGAATAGTTTAGACACTAGGATGAGATAAGGAACACACCCATTCT
AAAGAAATCACATTAGGATTCTCTTTTTAAGCTGTTCCTTAAAACACTAGAGTCTTAGAA
ATCTATTGGAGGCAGAAGCAGTCAAGGGTAGCCTAGGGTTAGGGTTAGGCTTAGGGTTAG
GGTTAGGGTACGGCTTAGGGTACTGTTTCGGGGAGGGGTTCAGGTACGGCGTAGGGTATG
GGTTAGGGTTAGGGTTAGGGTTAGTGTTAGGGTTAGGGCTCGGTTTAGGGTACGGGTTAG
GATTAGGGTACGTGTTAGGGTTAGGGTAGGGCTTAGGGTTAGGGTACGTGTTAGGGTTAG
GG
>gi:26|X53814.1|Blue Whale heavy satellite DNA
TAGTTATTAAACCTATCCCACTCTCTAGATACACCTTAGCACGTAAAGGAATATTATTTG
GGGGTCCAGACATGGAGAAGAGTTTAGACACTAGGATAAGATAAGGAACACACCCATTCT
AAAGAAATCACATTAGGATTCTCTTTTTAAGCTGTTCCTTAAAACTCTAGTGCTTAGGAA
ATCTATTGGAGGCAGAAGCAGTCAAGGGTAGCCTAGGGTTAGGGTTAGGCTTATGGTTAG
GGCTAGGGTACGGCTTAGGGTACGGATTCGGGGAGGGGTTCGGGTACGGCGTAGGGTATG
GGTTAGGGTTAGCGTTAGTGTTAGGGTTAGGGCTCGGTTTAGGGTACGGGTTAGGATTAG
GGTACGTGTTAGGGTTAGGGTAGGGGTTAGGGTTAGGGTACGCGTTAGGGTTAGGG
>gi:27|Z18633.1|B.physalus gene for large subunit rRNA
AACCAGTATTAGAGCACTGCCTGCCCGGTGACTAATCGTTAAACGGCCGCGGTATCCTGA
CCGTGCAAAGGTAGCATAATCACTTGTTCTCTAATTAGGGACTTGTATGAATGGCCACAC
GAGGGTTTTACTGTCTCTTACTTTTAATCAGTGAAATTGACCTCTCCGTGAAGAGGCGGA
GATAACAAAATAAGACGAGAAGACCCTATGGAGCTTCAATTAATCAACCCAAAAACCATA
ACCTTAAACCACCAAGGGATAACAAAACCTTATATGGGCTGACAATTTCGGTTGGGGTGA
CCTCGGAGTACAAAAAACCCTCCGAGTGATTAAAACTTAGGCCCACTAGCCAAAGTACAA
TATCACTTATTGATCCAATCCTTTGATCAACGGAACAAGTTACCCTAGGGATAACAGCGC
AATCCTATTCTAGAGTCCATATCGACAATAGGGTTTACGACCTCGATGTTGGATCAGGAC
ATCCTAATGGTGCAGCTGCTATTAAGGGTTCGTTTGTT
That's it,


Pierre

31 August 2012

The 500th post: Generating a pipeline of analysis (Makefile) for NGS with xslt.

This is my 500th post ! Happy birthday my blog ! (I know, I know, I have been unfaithful since the beginning of the year).

In my previous post, I've shown how a the .INTERMEDIATE keyword can be used in a Makefile to reduce the number of temporary files generated in a classical NGS pipeline. 'Make' also has an option '-j' to specify the number of threads: some independent tasks (such as aligning two distinct samples) can be processed in parallel .

At this point, I should cite this tweet from Sean Davis:



... and this even-more-recent tweet from Justin H. Johnson:




Here, I will show how XSLT can be used to generate a whole Makefile to process the output of Illumina/Casava after the FASTQs files have been generated using configureBclToFastq.pl using --fastq-cluster-count.

A Xml descriptor generated by Illumina Casava looks like this:
<DemultiplexConfig>
  <Software/>
  <FlowcellInfo ID="C0KTCACXX" Operator="JohnDoe" Recipe="" Desc="">
    <Lane Number="1">
      <Sample ProjectId="WGS" Control="N" Index="CTTGTA" SampleId="Father" Desc="1.5pM" Ref="" />
      <Sample ProjectId="WGS" Control="N" Index="GCCAAT" SampleId="Mother" Desc="1.5pM" Ref="" />
      <Sample ProjectId="WGS" Control="N" Index="TGACCA" SampleId="Daughter" Desc="1.5pM" Ref="" />
      <Sample ProjectId="Undetermined_indices" Control="N" Index="Undetermined" SampleId="lane1" Desc="Clusters with unmatched barcodes for lane 1" Ref="unknown" />
    </Lane>
    <Lane Number="2">
      <Sample ProjectId="WGS" Control="N" Index="CTTGTA" SampleId="Father" Desc="1.5pM" Ref="" />
      <Sample ProjectId="WGS" Control="N" Index="GCCAAT" SampleId="Mother" Desc="1.5pM" Ref="" />
      <Sample ProjectId="WGS" Control="N" Index="TGACCA" SampleId="Daughter" Desc="1.5pM" Ref="" />
      <Sample ProjectId="Undetermined_indices" Control="N" Index="Undetermined" SampleId="lane2" Desc="Clusters with unmatched barcodes for lane 2" Ref="unknown" />
    </Lane>
  </FlowcellInfo>
</DemultiplexConfig>
This document contains the required informations to find the fastqs and to generate a Makefile. I wrote the following XSLT stylesheet:


This stylesheet transforms the DemultiplexConfig.xml file to the following Makefile:

#
#files must have been generated using CASAVA/configureBclToFastq.pl with --fastq-cluster-count 0
#path to bwa
override BWA=/path/to/bwa-0.6.1/bwa
BWASAMPEFLAG= -a 600
BWAALNFLAG= -t 4
#path to samtools
SAMTOOLS=/path/to/samtools
#path to BEDTOOLS
BEDTOOLS=/path/to/bedtools/bin
#output directory
OUTDIR=2012831.OUTPUT.idp0
INPUTDIR=.
REF=/path/to/human_g1k_v37.fasta
JAVA=java

PICARDLIB=/path/to/picard-tools
TABIX=/path/to/tabix
GATK=/path/to/GenomeAnalysisTK.jar
GATKFLAGS=
VCFDBSNP=/path/to/ncbi/snp/00-All.vcf.gz

#rule how to create a bai from a bam:
%.bam.bai: %.bam
 #indexing BAM $<
 $(SAMTOOLS) index $<

#creates the default output directory
$(OUTDIR):
 mkdir -p $@

#indexes the reference for bwa
$(REF).bwt $(REF).ann $(REF).amb $(REF).pac $(REF).sa : $(REF)
 $(BWA) index -a bwtsw $<

#indexes the reference for samtools
$(REF).fai: $(REF)
 $(SAMTOOLS) faidx $<

#creates sequence dict for picard
$(REF).dict: $(REF)
 $(JAVA) -jar $(PICARDLIB)/CreateSequenceDictionary.jar REFERENCE=$(REF) OUTPUT=$<



#intermediate see http://plindenbaum.blogspot.fr/2012/08/next-generation-sequencing-gnu-make-and.html
.INTERMEDIATE :  \
 $(OUTDIR)/Project_WGS/Sample_Father/Father.bam.bai \
 $(OUTDIR)/Project_WGS/Sample_Father/Father.bam \
 $(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001.bam \
 $(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002.bam \
 $(OUTDIR)/Project_WGS/Sample_Mother/Mother.bam.bai \
 $(OUTDIR)/Project_WGS/Sample_Mother/Mother.bam \
 $(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001.bam \
 $(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002.bam \
 $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.bam.bai \
 $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.bam \
 $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001.bam \
 $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002.bam


#main target : generate the VCFs
all : vcfs


vcfs:  \
 $(OUTDIR)/Project_WGS/Sample_Father/Father.vcf.gz \
 $(OUTDIR)/Project_WGS/Sample_Mother/Mother.vcf.gz \
 $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.vcf.gz


$(OUTDIR)/Project_WGS/Sample_Father/Father.realigned.bam : $(OUTDIR)/Project_WGS/Sample_Father/Father.bam $(OUTDIR)/Project_WGS/Sample_Father/Father.bam.bai
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T RealignerTargetCreator \
  -R $(REF) \
  -I $< \
  -o $<.intervals \
  --known $(VCFDBSNP)
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T IndelRealigner \
  -R $(REF) \
  -I $< \
  -o $@ \
  -targetIntervals $<.intervals \
  --knownAlleles $(VCFDBSNP)
 rm -f $<.intervals


#creates a BAM for 'Father' by merging the other lanes
$(OUTDIR)/Project_WGS/Sample_Father/Father.bam: $(OUTDIR)  \
  $(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001.bam \
  $(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002.bam
 #merge BAMS
 $(JAVA) -jar $(PICARDLIB)/MergeSamFiles.jar SORT_ORDER=coordinate OUTPUT=$@ \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001.bam \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002.bam
 #fix mate information
 $(JAVA) -jar $(PICARDLIB)/FixMateInformation.jar \
  INPUT=$@
 #validate
 $(JAVA) -jar $(PICARDLIB)/ValidateSamFile.jar \
  VALIDATE_INDEX=true \
  I=$@


#Creates a BAM for 'Project_WGS/Sample_Father/Father_CTTGTA_L001'

$(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Father
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp22528_1.sai $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp22528_2.sai $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idp22528 CN:MyLaboratory LB:Father PG:bwa SM:Father PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idp22528.sam \
  $(OUTDIR)/idp22528_1.sai \
  $(OUTDIR)/idp22528_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L001_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idp22528_1.sai $(OUTDIR)/idp22528_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idp22528_unsorted.bam $(OUTDIR)/idp22528.sam
 #remove the sam
 rm -f $(OUTDIR)/idp22528.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idp22528_unsorted.bam $(OUTDIR)/idp22528_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idp22528_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idp22528_sorted.bam \
  OUTPUT=$(OUTDIR)/idp22528_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Father/Father.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idp22528_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idp22528_dup.bam $@



$(OUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Father
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp31344_1.sai $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp31344_2.sai $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idp31344 CN:MyLaboratory LB:Father PG:bwa SM:Father PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idp31344.sam \
  $(OUTDIR)/idp31344_1.sai \
  $(OUTDIR)/idp31344_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Father/Father_CTTGTA_L002_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idp31344_1.sai $(OUTDIR)/idp31344_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idp31344_unsorted.bam $(OUTDIR)/idp31344.sam
 #remove the sam
 rm -f $(OUTDIR)/idp31344.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idp31344_unsorted.bam $(OUTDIR)/idp31344_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idp31344_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idp31344_sorted.bam \
  OUTPUT=$(OUTDIR)/idp31344_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Father/Father.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idp31344_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idp31344_dup.bam $@




#creates a VCF for 'Father'
$(OUTDIR)/Project_WGS/Sample_Father/Father.vcf.gz : $(OUTDIR)/Project_WGS/Sample_Father/Father.realigned.bam
 $(JAVA) -jar $(GATK) \
  -T UnifiedGenotyper \
  -glm BOTH \
  -baq OFF \
  -R $(REF) \
  -I $< \
  --dbsnp $(VCFDBSNP) \
  -o $(OUTDIR)/Project_WGS/Sample_Father/Father.vcf
 $(TABIX)/bgzip $(OUTDIR)/Project_WGS/Sample_Father/Father.vcf
 $(TABIX)/tabix -f -p vcf $(OUTDIR)/Project_WGS/Sample_Father/Father.vcf.gz
$(OUTDIR)/Project_WGS/Sample_Mother/Mother.realigned.bam : $(OUTDIR)/Project_WGS/Sample_Mother/Mother.bam $(OUTDIR)/Project_WGS/Sample_Mother/Mother.bam.bai
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T RealignerTargetCreator \
  -R $(REF) \
  -I $< \
  -o $<.intervals \
  --known $(VCFDBSNP)
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T IndelRealigner \
  -R $(REF) \
  -I $< \
  -o $@ \
  -targetIntervals $<.intervals \
  --knownAlleles $(VCFDBSNP)
 rm -f $<.intervals


#creates a BAM for 'Mother' by merging the other lanes
$(OUTDIR)/Project_WGS/Sample_Mother/Mother.bam: $(OUTDIR)  \
  $(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001.bam \
  $(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002.bam
 #merge BAMS
 $(JAVA) -jar $(PICARDLIB)/MergeSamFiles.jar SORT_ORDER=coordinate OUTPUT=$@ \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001.bam \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002.bam
 #fix mate information
 $(JAVA) -jar $(PICARDLIB)/FixMateInformation.jar \
  INPUT=$@
 #validate
 $(JAVA) -jar $(PICARDLIB)/ValidateSamFile.jar \
  VALIDATE_INDEX=true \
  I=$@


#Creates a BAM for 'Project_WGS/Sample_Mother/Mother_GCCAAT_L001'

$(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Mother
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp15352_1.sai $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp15352_2.sai $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idp15352 CN:MyLaboratory LB:Mother PG:bwa SM:Mother PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idp15352.sam \
  $(OUTDIR)/idp15352_1.sai \
  $(OUTDIR)/idp15352_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L001_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idp15352_1.sai $(OUTDIR)/idp15352_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idp15352_unsorted.bam $(OUTDIR)/idp15352.sam
 #remove the sam
 rm -f $(OUTDIR)/idp15352.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idp15352_unsorted.bam $(OUTDIR)/idp15352_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idp15352_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idp15352_sorted.bam \
  OUTPUT=$(OUTDIR)/idp15352_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Mother/Mother.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idp15352_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idp15352_dup.bam $@



$(OUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Mother
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idm22432_1.sai $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idm22432_2.sai $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idm22432 CN:MyLaboratory LB:Mother PG:bwa SM:Mother PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idm22432.sam \
  $(OUTDIR)/idm22432_1.sai \
  $(OUTDIR)/idm22432_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Mother/Mother_GCCAAT_L002_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idm22432_1.sai $(OUTDIR)/idm22432_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idm22432_unsorted.bam $(OUTDIR)/idm22432.sam
 #remove the sam
 rm -f $(OUTDIR)/idm22432.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idm22432_unsorted.bam $(OUTDIR)/idm22432_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idm22432_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idm22432_sorted.bam \
  OUTPUT=$(OUTDIR)/idm22432_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Mother/Mother.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idm22432_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idm22432_dup.bam $@




#creates a VCF for 'Mother'
$(OUTDIR)/Project_WGS/Sample_Mother/Mother.vcf.gz : $(OUTDIR)/Project_WGS/Sample_Mother/Mother.realigned.bam
 $(JAVA) -jar $(GATK) \
  -T UnifiedGenotyper \
  -glm BOTH \
  -baq OFF \
  -R $(REF) \
  -I $< \
  --dbsnp $(VCFDBSNP) \
  -o $(OUTDIR)/Project_WGS/Sample_Mother/Mother.vcf
 $(TABIX)/bgzip $(OUTDIR)/Project_WGS/Sample_Mother/Mother.vcf
 $(TABIX)/tabix -f -p vcf $(OUTDIR)/Project_WGS/Sample_Mother/Mother.vcf.gz
$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.realigned.bam : $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.bam $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.bam.bai
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T RealignerTargetCreator \
  -R $(REF) \
  -I $< \
  -o $<.intervals \
  --known $(VCFDBSNP)
 $(JAVA) -jar $(GATK) $(GATKFLAGS) \
  -T IndelRealigner \
  -R $(REF) \
  -I $< \
  -o $@ \
  -targetIntervals $<.intervals \
  --knownAlleles $(VCFDBSNP)
 rm -f $<.intervals


#creates a BAM for 'Daughter' by merging the other lanes
$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.bam: $(OUTDIR)  \
  $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001.bam \
  $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002.bam
 #merge BAMS
 $(JAVA) -jar $(PICARDLIB)/MergeSamFiles.jar SORT_ORDER=coordinate OUTPUT=$@ \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001.bam \
  INPUT=$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002.bam
 #fix mate information
 $(JAVA) -jar $(PICARDLIB)/FixMateInformation.jar \
  INPUT=$@
 #validate
 $(JAVA) -jar $(PICARDLIB)/ValidateSamFile.jar \
  VALIDATE_INDEX=true \
  I=$@


#Creates a BAM for 'Project_WGS/Sample_Daughter/Daughter_TGACCA_L001'

$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Daughter
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp5176_1.sai $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp5176_2.sai $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idp5176 CN:MyLaboratory LB:Daughter PG:bwa SM:Daughter PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idp5176.sam \
  $(OUTDIR)/idp5176_1.sai \
  $(OUTDIR)/idp5176_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L001_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idp5176_1.sai $(OUTDIR)/idp5176_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idp5176_unsorted.bam $(OUTDIR)/idp5176.sam
 #remove the sam
 rm -f $(OUTDIR)/idp5176.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idp5176_unsorted.bam $(OUTDIR)/idp5176_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idp5176_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idp5176_sorted.bam \
  OUTPUT=$(OUTDIR)/idp5176_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idp5176_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idp5176_dup.bam $@



$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002.bam: $(REF).bwt \
 $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R1_001.fastq.gz \
 $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R2_001.fastq.gz
 #create directory
 mkdir -p $(OUTDIR)/Project_WGS/Sample_Daughter
 #align fastq1 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp211744_1.sai $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R1_001.fastq.gz
 #align fastq2 with BWA
 $(BWA) aln $(BWAALNFLAG) $(REF) -f $(OUTDIR)/idp211744_2.sai $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R2_001.fastq.gz
 #sampe those files
 $(BWA) sampe $(BWASAMPEFLAG)  -r '@RG ID:idp211744 CN:MyLaboratory LB:Daughter PG:bwa SM:Daughter PL:ILLUMINA' $(REF) \
  -f $(OUTDIR)/idp211744.sam \
  $(OUTDIR)/idp211744_1.sai \
  $(OUTDIR)/idp211744_2.sai \
  $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R1_001.fastq.gz \
  $(INPUTDIR)/Project_WGS/Sample_Daughter/Daughter_TGACCA_L002_R2_001.fastq.gz
 #remove sai
 rm -f $(OUTDIR)/idp211744_1.sai $(OUTDIR)/idp211744_2.sai
 #convert sam to bam
 $(SAMTOOLS) view -b -T $(REF) -o $(OUTDIR)/idp211744_unsorted.bam $(OUTDIR)/idp211744.sam
 #remove the sam
 rm -f $(OUTDIR)/idp211744.sam
 #sort this bam
 $(SAMTOOLS) sort  $(OUTDIR)/idp211744_unsorted.bam $(OUTDIR)/idp211744_sorted
 #remove the unsorted bam
 rm -f $(OUTDIR)/idp211744_unsorted.bam
 #mark duplicates
 $(JAVA) -jar $(PICARDLIB)/MarkDuplicates.jar VALIDATION_STRINGENCY=LENIENT ASSUME_SORTED=true TMP_DIR=$(OUTDIR) \
  INPUT=$(OUTDIR)/idp211744_sorted.bam \
  OUTPUT=$(OUTDIR)/idp211744_dup.bam \
  METRICS_FILE=$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.metrics
 #remove before duplicate
 rm -f $(OUTDIR)/idp211744_sorted.bam
 #move to final bam
 mv $(OUTDIR)/idp211744_dup.bam $@




#creates a VCF for 'Daughter'
$(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.vcf.gz : $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.realigned.bam
 $(JAVA) -jar $(GATK) \
  -T UnifiedGenotyper \
  -glm BOTH \
  -baq OFF \
  -R $(REF) \
  -I $< \
  --dbsnp $(VCFDBSNP) \
  -o $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.vcf
 $(TABIX)/bgzip $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.vcf
 $(TABIX)/tabix -f -p vcf $(OUTDIR)/Project_WGS/Sample_Daughter/Daughter.vcf.gz
(Of course, this workflow would satisfy my needs, may be not yours)

Et voila !


That's it !


Pierre