EricChu:LabNotesMDA/2014-4-16: Difference between revisions
>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. | 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. | 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) | ||
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=",")) | > y8<-c(scan("PGP1#21CoREChamber8.csv", sep=",")) | ||
Read 50000 items | Read 50000 items | ||
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 ... | ||
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) | ||
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: | > 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 | ||
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. | |||
[[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.
- 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.