EricChu:LabNotesMDA/2014-4-16

From ZhangLabWiki
Jump to navigation Jump to search

Fragment Calling Protocol[edit]

  • R code smooth.discrete in the mhsmm package

1. Starting from the BAM files, convert to SAM from BAM

samtools view -o s_9_Indx22_unique.sam s_9_Indx22_unique.bam

2. Python counting reads with 50k varbin

python varbin.50k.sam.py s_9_Indx22_unique.sam s_9_Indx22.varbin.50k.txt s_9_Indx22.varbin.50k.stats.txt

3. Open the varbin.50k file in excel. It has 50,000 rows. Column 4 is the read count of that bin.

4. Copy column 4 to a new sheet

5. Use a function to turn (read counts>20 to a value of 2 and otherwise =1). “2” represents the presence of a DNA strand.

=IF(B1>=20,2,1)		

6. Only save the converted values in CSV form.

7. Repeat this for each chamber for analysis.

8. Load the mshmm package to R

9. Input the CSV file in R as a vector

> y8<-c(scan("PGP1#21CoREChamber8.csv", sep=","))
Read 50000 items

10. Make sure the size was right

> str(y8)
num [1:50000] 2 2 1 1 1 1 1 1 1 1 ...

11. Set parameters

> initial=c(1,0)
> P=matrix(c(0.99,0.01,0.01,0.99),nrow=2)   #transition matrix
> B=P   #emission matrix
# The matrices will eventually be trained and then reflect the actual probability.
# as long as the transition matrix was not too far from the reality (>90% chance of no change)

12. Run all vectors through smooth.discrete (HMM) and plot

> obj8=smooth.discrete(y8,init=initial,trans=P,parms.emission=B)
> plot (y8, ylim=c(0.8,2))
> addStates(obj8$s)
# the addStates function will add a colored bar at the bottom of the plot
> plot (y30[2800:3000], ylim=c(0.8,2)) #plot a specific region
> plot(y30, ylim=c(0.8,2)) # or plot the entire genome
> addStates(obj30$s[2800:3000]) # add the visualization of this specific region
> addStates(obj30$s) # or add to the entire genome

13. The state is stored in the vector “obj30$s”. It is the predicted states of 1s and 2s in a vector form.

14. Output the state file to obj30.csv

> write.table(t(obj30$s), "obj30.csv", row.names=FALSE, col.names=FALSE, sep=",")

15. Finding the transition and emission matrices after training

> summary (obj30)
init: 
 0 1 

transition:
      [,1]  [,2]
[1,] 0.969 0.031
[2,] 0.182 0.818

emission:
$pmf
           [,1]       [,2]
[1,] 0.97018569 0.02981431
[2,] 0.04952785 0.95047215


Compare Region 165Mbp-169Mbp of Chromosome 4 from the chambers 10,14,22,24 of PGP1#21 CoRE[edit]

  • Figure below showing region 160M-170Mbp of Chromosome 4. There was no leakage between neighboring chambers.

File:PGP1 21CoREAllChambersofChr4.jpg

  • Figure focusing on the region 164M-170Mbp of Chromosome 4 in 4 chambers 10,14,22,24 (top to bottom).

File:PGP1 12CoRE4ChambersofChr4.jpg

  • Superimpose the Varbin.50k of the same region to the R prediction map. We can find a pretty nice match of the expected fragments being called.

File:PGP1 21CoREChamber10.jpg

File:PGP1 21CoREChamber14.jpg

File:PGP1 21CoREChamber22.jpg

File:PGP1 21CoREChamber24.jpg

  • Directly compare the prediction of these 4 chambers. Clear overlaps in these regions were found.

File:PGP1 21CoREVarbin Prediction4chambers.jpg

  • Stacking up chambers 9-24. Orange represents the region with DNA coverage.
  • By combining these vectors, pulling each bin with "2" and calling this region covered. Then taking the average of all these vectors.

File:PGP1 21CoREVarbin PredictionWholeCell.jpg