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

From ZhangLabWiki
Jump to navigation Jump to search
>Ericchu
(Created page with "==Fragment Calling Protocol== * 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. 2...")
 
>Ericchu
No edit summary
 
(6 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
2. Copy column 4 to a new sheet
samtools view -o s_9_Indx22_unique.sam s_9_Indx22_unique.bam
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.  
 
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)
  =IF(B1>=20,2,1)
4. Only save the converted values in CSV form.
 
5. Repeat this for each chamber for analysis.
6. Only save the converted values in CSV form.
6. Load the mshmm package to R
 
7. Input the CSV file in R as a vector
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=","))
  > 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 20: Line 35:
  # The matrices will eventually be trained and then reflect the actual probability.
  # 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)
  # 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 26: 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