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

From ZhangLabWiki
Jump to navigation Jump to search
>Ericchu
No edit summary
>Ericchu
No edit summary
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 40: Line 47:
  > 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.
 




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

Revision as of 21:18, 17 April 2014

Fragment Calling Protocol

  • 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: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

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