jimmy大神的教程对我还是有点难度,主要是直接操作xshell水平不太够
使用xshell+winSCP操作
一、数据下载
wget sra数据网址
#拆分
fastq-dump -- gzip -- split - files SRR
二、质控
#!/bin/sh
fastp -i test_R1.fastq -I test_R2.fastq -o outtest_R1.fq.gz -O outtest_R2.fq.gz
三、比对
#!/bin/sh
bowtie2-build -f genomic.fna.gz --threads 24 pig(物种名)
bowtie2 -p 10 -x ../bidui/pig -1 outtest_R1.fq.gz -2 /outtest_R2.fq.gz -S test.sam
samtools view -F 4 -Su test.sam | samtools sort -T 11111 -o test_sort.bam
samtools index test_sort.bam
#(sort.bam为排序后的bam)
#查看比对率
#!/bin/sh
qualimap --java-mem-size=50G bamqc -bam test.bam -outdir ../ -outformat PDF:HTML
四、去除线粒体、PCR重复等
#!/bin/sh
gatk MarkDuplicates \
-I test.bam -O test.picard.redup.bam --REMOVE_SEQUENCING_DUPLICATES true -M test.log
samtools index test.picard.redup.bam
samtools flagstat test.picard.redup.bam > test.rmdup.stat
samtools view -h -f 2 -q 30 test.picard.redup.bam |grep -v chrM|samtools sort -O bam -@ 5 -o - >test.last.bam
samtools index test.last.bam
samtools flagstat test.last.bam > test.last.stat
bedtools bamtobed -i test.last.bam > test.bed
五、peaks calling
#!/bin/sh
ls *.bed | while read id ;do(/nfs3/wanghlab/anaconda3/envs/chip/bin/macs2 callpeak -t $id -g mm --nomodel --shift -100 --extsize 200 - n ${id%%.*} --outdir ../);done
后续步骤晚些再写