EricChu:LabNotesMDA/2014-4-16: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Ericchu
No edit summary
>Ericchu
No edit summary
 
(5 intermediate revisions by the same user not shown)
Line 1: Line 1:
==Fragment Calling Protocol==
==Fragment Calling Protocol==
* R code smooth.discrete in the mhsmm package
* R code smooth.discrete in the mhsmm package
1. Open the varbin.50k file in excel. It has 50,000 rows. Column 4 is the read count of that bin.
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. Copy column 4 to a new sheet
2. Python counting reads with 50k varbin


3. Use a function to turn (read counts>20 to a value of 2 and otherwise =1). “2” represents the presence of a DNA strand.  
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)
  =IF(B1>=20,2,1)


4. Only save the converted values in CSV form.
6. Only save the converted values in CSV form.


5. Repeat this for each chamber for analysis.
7. Repeat this for each chamber for analysis.


6. Load the mshmm package to R
8. Load the mshmm package to R


7. Input the CSV file in R as a vector
9. Input the CSV file in R as a vector
  > y8<-c(scan("PGP1#21CoREChamber8.csv", sep=","))
  > y8<-c(scan("PGP1#21CoREChamber8.csv", sep=","))
  Read 50000 items
  Read 50000 items


8. Make sure the size was right
10. Make sure the size was right
  > str(y8)
  > str(y8)
  num [1:50000] 2 2 1 1 1 1 1 1 1 1 ...
  num [1:50000] 2 2 1 1 1 1 1 1 1 1 ...


9. Set parameters
11. Set parameters
  > initial=c(1,0)
  > initial=c(1,0)
  > P=matrix(c(0.99,0.01,0.01,0.99),nrow=2)  #transition matrix
  > P=matrix(c(0.99,0.01,0.01,0.99),nrow=2)  #transition matrix
Line 29: Line 36:
  # as long as the transition matrix was not too far from the reality (>90% chance of no change)
  # as long as the transition matrix was not too far from the reality (>90% chance of no change)


10. Run all vectors through smooth.discrete (HMM) and plot
12. Run all vectors through smooth.discrete (HMM) and plot
  > obj8=smooth.discrete(y8,init=initial,trans=P,parms.emission=B)
  > obj8=smooth.discrete(y8,init=initial,trans=P,parms.emission=B)
  > plot (y8, ylim=c(0.8,2))
  > plot (y8, ylim=c(0.8,2))
Line 35: Line 42:
  # the addStates function will add a colored bar at the bottom of the plot
  # the addStates function will add a colored bar at the bottom of the plot


  > plot (y30[2800:300], ylim=c(0.8,2)) #plot a specific region
  > 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
  > 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[2800:3000]) # add the visualization of this specific region
  > addStates(obj30$s) # or add to the entire genome
  > addStates(obj30$s) # or add to the entire genome


11. The state is stored in the vector “obj30$s”. It is the predicted states of 1s and 2s in a vector form.
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==
* Figure below showing region 160M-170Mbp of Chromosome 4. There was no leakage between neighboring chambers.
 
[[File:PGP1_21CoREAllChambersofChr4.jpg| 1000px]]
 
* Figure focusing on the region 164M-170Mbp of Chromosome 4 in 4 chambers 10,14,22,24 (top to bottom).
 
[[File:PGP1_12CoRE4ChambersofChr4.jpg| 1000px]]
 
* 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| 500px]]
 
[[File:PGP1_21CoREChamber14.jpg| 500px]]
 
[[File:PGP1_21CoREChamber22.jpg| 500px]]
 
[[File:PGP1_21CoREChamber24.jpg| 500px]]
 
* Directly compare the prediction of these 4 chambers. Clear overlaps in these regions were found.


[[File:PGP1_21CoREVarbin_Prediction4chambers.jpg| 500px]]


* 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.


* Compare Region 165Mbp-169Mbp of Chromosome 4 from the chambers 10,14,22,24 of PGP1#21 CoRE
[[File:PGP1_21CoREVarbin_PredictionWholeCell.jpg| 500px]]

Latest revision as of 00:30, 24 April 2014

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