Our own data alignment and analysis: Difference between revisions

From ZhangLabWiki
Jump to navigation Jump to search
>Shicheng
(Created page with " === fastq merge by indx === 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 & d...")
 
>Shicheng
 
(17 intermediate revisions by the same user not shown)
Line 1: Line 1:


=== Transfer data from Genome-miner to TSCC===
* md5 to check the integrity
cd /home/shg047/oasis/Alice/mouse/bam/iPScnt/bam
for i in `ls Indx*bam`
do
samtools fastq $i > $i.fastq &
done
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 ===
=== 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}
  for i in {01,02,04,05,06,07,09,10,11,12}
Line 8: Line 53:
  cat s_*_2_Indx$i.txt.gz > Indx$i.read2.fq.gz &
  cat s_*_2_Indx$i.txt.gz > Indx$i.read2.fq.gz &
  done
  done
=== RD90 and UP80 ===
[[RD90 and L80bp Regions]]
* 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 ===
=== 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 ===


=== bismark alignment ===
=== 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

Latest revision as of 08:15, 26 January 2017


Transfer data from Genome-miner to TSCC[edit]

  • md5 to check the integrity
cd /home/shg047/oasis/Alice/mouse/bam/iPScnt/bam
for i in `ls Indx*bam`
do
samtools fastq $i > $i.fastq &
done
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[edit]

  • 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[edit]

RD90 and L80bp Regions

  • 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[edit]

bismark alignment[edit]

scp shg047@genome-miner.ucsd.edu:/media/Home_Raid1/zhl002/NAS2/LTS33T/mouse_WGBS/BAMfiles/*  /home/shg047/oasis/mouse/alice/bam

RD>10 Regions[edit]

Random Sampling Haploinfo by Bootstrap[edit]

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