Our own data alignment and analysis
Jump to navigation
Jump to search
Transfer data from Genome-miner to TSCC
- md5 to check the integrity
cd /home/shg047/oasis/mouse/alice/raw
for i in `ls *txt.gz` do md5sum $i > $i.md5 & done
scp shg047@genome-miner.ucsd.edu:/media/NAS3_volume1/SeqStore2016/130104_SN1001/* ./
for i in `ls *txt.gz` do md5sum $i > ~/$i.md5 & done
fastq merge by indx
- Trim fastq separately
my @file=glob("*gz"); my %file; foreach my $file(@file){ my @line=split/\.|_/,$file; push @{$file{$line[1].$line[3]}},$file; } foreach my $id(sort keys %file){ foreach my $sam(@{$file{$id}}){ print "$sam\t"; } print "\n"; } perl compile.pl > FastMatch.txt perl ~/bin/trimGolrePBS.pl FastMatch.txt submit
- bismark alignment separately
for i in {01,02,04,05,06,07,09,10,11,12} do cat s_*_1_Indx$i.txt.gz > Indx$i.read1.fq.gz & cat s_*_2_Indx$i.txt.gz > Indx$i.read2.fq.gz & done
RD90 and UP80
- Merge Bams from Different samples
samtools merge -o output.bam $1 $2 $3
- Split big bam to small bam by chrosome
for i in {1..19} X Y M do samtools view -b -q 20 Mouse.MHB.Alice.MergeBam.sort.bam chr$i > Mouse.MHB.Alice.chr$i.bam & done
- Obtain depth for each postion
for i in `ls *bam` do bedtools genomecov -ibam $i -bg > merge.$i.bedcov done
- Obtain D>=30 and Len>80 regions
for i in {1..19} X Y M do cd /home/shg047/oasis/mouse/RD awk '$4>29 { print $1"\t"$2"\t"$3}' merge.chr$i.bam.gencov.bed | bedtools merge -d 10 -i - > merge.chr$i.bamRD10.bed awk '$3-$2>80 {print $1"\t"$2"\t"$3"\t"$3-$2+1}' merge.chr$i.bamRD10.bed > merge.chr$i.bam.RD10_80up.bed done
for i in chr10:79857968-79858049 chr10:79858088-79858205 chr10:79858211-79858298 chr10:79858319-79858347 chr6:122658641-122658674 chr6:122659252-122659297 do perl ~/bin/hapinfo2LDR2.pl $i.tissue.R2.txt $i < tissue.hapinfo.txt perl ~/bin/hapinfo2LDR2.pl $i.ips.R2.txt $i < ipsNT.hapinfo.txt perl ~/bin/hapinfo2LDR2.pl $i.esc.R2.txt $i < ESC.hapinfo.txt done
- LD heatmap plot
file<-list.files(pattern="*.rsq$") for(i in file){ print(i) M<-read.table(i,head=T,row.names=1,as.is=T) library("grDevices") col=colorRampPalette(c("white", "red"))(10) M[lower.tri(M)] <- NA pdf(paste(i,"pdf",sep=".")) image(data.matrix(M),col = col,frame=F,xaxt="n",yaxt="n") dev.off() }
Trim_galore
bismark alignment
scp shg047@genome-miner.ucsd.edu:/media/Home_Raid1/zhl002/NAS2/LTS33T/mouse_WGBS/BAMfiles/* /home/shg047/oasis/mouse/alice/bam
RD>10 Regions
Random Sampling Haploinfo by Bootstrap
for i in {1..100} do perl ~/bin/randomSampleFromHaploInfo.pl ipsNT.hapinfo.txt > ipsNT.hapinfo.txt.$i perl ~/bin/randomSampleFromHaploInfo.pl ESC.hapinfo.txt.$i perl ~/bin/randomSampleFromHaploInfo.pl tissue.hapinfo.txt.$i echo $i done