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
Line 2: Line 2:
* 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. Open the varbin.50k file in excel. It has 50,000 rows. Column 4 is the read count of that bin.
2. Copy column 4 to a new sheet
2. Copy column 4 to a new sheet
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.  
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.  
  =IF(B1>=20,2,1)
  =IF(B1>=20,2,1)
4. Only save the converted values in CSV form.
4. Only save the converted values in CSV form.
5. Repeat this for each chamber for analysis.
5. Repeat this for each chamber for analysis.
6. Load the mshmm package to R
6. Load the mshmm package to R
7. Input the CSV file in R as a vector
7. 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
8. 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
9. Set parameters
  > initial=c(1,0)
  > initial=c(1,0)
Line 20: Line 28:
  # 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
10. 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)
Line 30: Line 39:
  > 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.
11. The state is stored in the vector “obj30$s”. It is the predicted states of 1s and 2s in a vector form.



Revision as of 21:15, 17 April 2014

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. Copy column 4 to a new sheet

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.

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

4. Only save the converted values in CSV form.

5. Repeat this for each chamber for analysis.

6. Load the mshmm package to R

7. Input the CSV file in R as a vector

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

8. Make sure the size was right

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

9. 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)

10. 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:300], 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

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


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