EricChu:LabNotesMDA/2014-4-16: Difference between revisions
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