AlanFung:LabNotes/ASE/2010/2012-7-18

From ZhangLabWiki
Jump to navigation Jump to search

Capturing Protocol[edit]

[[1]]

Analysis of Gimelbrant Lab Samples[edit]

  • Data Stored in seqStore/120712_HL125, transfer my data to Gimelbrant/ase_06262012/tmp directory
  • Transfer to my own directory
  • Create folder gimelbrant/ase_06262012
  • Make a directory for each samples (bash mkdir.txt)

Combine files and transfer to individual folders[edit]

  • Run combine.txt to combine data and separate them into individual folders
  • Run move.txt to transfer files into individual folders

AseBowtieSam[edit]

  • Run aseBowtieSam.pl
  • chmod +x../aseBowtie.pl to make is executable
  • make sure the directories are up to date by running each line manually if not up to date edit the script
my $mapper = "/home/kunzhang/softwares/bowtie-0.11.3/bowtie";
* change to my $mapper = "/home/kunzhang/softwares/bowtie-0.12.8/bowtie";
my $samtool = "/home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/samtools";
my $readTrimmer="/home/kunzhang/Exome/Solexa/scripts/trimTerminalMismatches.pl";
*Cannot find this one
*change to my $readTrimmer="/home/kunzhang/Exome/Data/scripts/trimTerminalMismatches.pl";
my $bowtie2sam = "/home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl";
my $snpFilter = "/home/kunzhang/MIP/pileup2basecalls.pl";
my $in_file = $ARGV[0];
my $template = "/home/kunzhang/HsGenome/hg18.masked/hg18.masked";
my $fai_file = "/home/kunzhang/HsGenome/hg18.masked/HsGenome18.masked.fai";
  • Run the script with the .txt file
  • after editing the script, copy script to cDNA-DF2
  • run script
  • Return with error
/usr/local/bin/maq sol2sanger cDNA-DF2.txt cDNA-DF2.fastq
Can't exec "/usr/local/bin/maq": No such file or directory at aseBowtieSam.pl line 28.
Warning: Could not open read file "cDNA-DF2.fastq" for reading; skipping...
sh: cannot open cDNA-DF2.bowtie.out: No such file
Can't open cDNA-DF2.bowtie.out.trimmed: No such file or directory at /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl line 19.
Can't use an undefined value as an ARRAY reference at /home/kunzhang/softwares/samtools-0.1.7_x86_64-linux/bowtie2sam.pl line 48.
  • It appears that maq is uninstalled, Athurva will help me install the software

Shell to run script on all samples[edit]

*for f in gDNA-H14 ChIPK36me3-P10 ChIPK36me3-DF2 ChIPK27me3-P10 ChIPK27me3-DF2 cDNA-P10E3 cDNA-DF2 gDNA-GM gDNA-P10E3 gDNA-DF2 gDNA-H0 gDNA-GmgDNA:H0gDNA::1:1   
do
cd $f
../aseBowtieSam.pl $f.txt >$f.log
cd ..
done
  • Run is complete
  • What this script does? Looking at the output it generates
.bowtie.basecalls.txt
.pileup
.sorted.bam
  • Going back to Athurava's lecture notes
  • Pileup is the consensus sequence, it describes the base-pair information at each chromosomal position
Pileup File Columns
1. Chromosome
2. Position
3. reference base
4. Call
5. Consensus quality 

P(call column is incorrect) Coverage sum of quality score at that position

Analysis for each sample[edit]

Total reads[edit]

  • wc -l .txt >> .txt
  • go into the fastq .txt files and divide the total lines by 4
  • run it using script
  • copy to excel data to text column /4

Mapping Rate[edit]

  • /../home/kunzhang/softwares/samtools-latest/samtools flagstat cDNA-DF2.sorted.bam
  • line 4 is the mapped reads
  • sed -n 4p $f.log2 >> ../mapreadcount.txt
  • to print line 4 into text file

#SNPs called[edit]

  • wc -l $f.bowtie.basecalls.txt >> ../snpcount.txt
  • Basically count the lines in the basecalls.txt file

#Het called[edit]

  • Kun told me i can use the IF function in excel and count the heterozygous SNP
  • Used IF function for column C , =IF(LEFT(C1)=RIGHT(C1),0,1)
  • If the two base are the same returns 0
  • If the two bases are different returns 1
  • Get the total of the column and you get the total number of heterozygous SNP

Data Summary[edit]

Sample Total Number of Reads Mapping Reads Mapping Rate #SNPs called #Het % Het
gDNA-H7 7,214,167 4,946,303 68.56 28,280 2954 10.45
gDNA-H14 9,686,002 8,221,702 84.88 20,536 2502 12.18
ChIPK36me3-P10 8,811,298 7,223,409 81.98 25,452 2451 9.63
ChIPK36me3-DF2 9,243,720 7,649,171 82.75 26,507 2660 10.04
ChIPK27me3-P10 9,149,672 7,490,869 81.87 25,320 2475 9.77
ChIPK27me3-DF2 8,125,103 6,607,398 81.32 26,137 2687 10.28
cDNA-P10E3 6,012,037 5,207,362 86.62 19,096 2451 12.84
cDNA-DF2 4,306,791 3,756,778 87.23 18,854 2405 12.76
gDNA-GM 6,920,955 5,504,556 79.53 30,270 4188 13.84
gDNA-P10E3 7,540,086 6,126,998 81.26 30,887 4223 13.67
gDNA-DF2 6,836,212 5,328,454 77.94 30,249 4148 13.71
gDNA-H0 8,268,386 6,506,594 78.69 30,809 3590 11.65
gDNA-GmgDNA:H0gDNA::1:1    9,113,177 7,336,107 80.50 29,753 4065 13.66