转录组数据比对subjunc-7

## ----构建索引# 进入参考基因组目录cd $HOME/database/GRCh38.105

# subjunc构建索引，构建索引时间比较长大约40分钟左右，建议携程sh脚本提交后台运行# vim SubjuncIndex.shmkdir SubjuncIndexfa=Homo_sapiens.GRCh38.dna.primary_assembly.fafa_baseName=GRCh38.dnasubread-buildindex -o SubjuncIndex/${fa_baseName} ${fa}

# 运行nohup sh SubjuncIndex.sh >SubjuncIndex.sh.log &

## ----比对# 进入文件夹cd $HOME/project/Human-16-Asthma-Trans/Mapping/Subjunc

# vim subjunc.shindex=/home/t_rna/database/GRCh38.104/SubjuncIndex/GRCh38.dnainputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galoreoutdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Subjunc cat ../../data/cleandata/trim_galore/ID | while read iddo subjunc -T 5 -i ${index} -r ${inputdir}/${id}_1_val_1.fq.gz -R ${inputdir}/${id}_2_val_2.fq.gz -o ${outdir}/${id}.Subjunc.bam 1>${outdir}/${id}.Subjunc.log 2>&1 && samtools sort -@ 6 -o ${outdir}/${id}.Subjunc.sorted.bam ${outdir}/${id}.Subjunc.bamdone # 运行nohup sh subjunc.sh >subjunc.log &

# 单个样本samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam # 多个样本，vim flagstat.shls *.sorted.bam | while read iddo samtools flagstat -@ 3 ${id} > ${id/bam/flagstat}done# 质控multiqc -o ./ *.flagstat # 运行nohup sh flagstat.sh >flagstat.log &

##----view查看bam文件samtools view SRR1039510.Hisat_aln.sorted.bamsamtools view -H SRR1039510.Hisat_aln.sorted.bamsamtools view -h SRR1039510.Hisat_aln.sorted.bam ##----index对bam文件建索引samtools index SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sorted.bam.bai ##----flagstat统计比对结果samtools flagstat -@ 3 SRR1039510.Hisat_aln.sorted.bam ##----sort排序 sam转bam并排序samtools sort -@ 3 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam ##----depth统计测序深度# 得到的结果中，一共有3列以指标分隔符分隔的数据，第一列为染色体名称，第二列为位点，第三列为覆盖深度samtools depth SRR1039510.Hisat_aln.sorted.bam >SRR1039510.Hisat_aln.sorted.bam.depth.txt ##----计算某一个基因的测序深度# 找到基因的坐标zless -S Homo_sapiens.GRCh38.95.gff3.gz |awk '{if($3=="gene")print}' |grep 'ID=gene:ENSG00000186092' |awk '{print $1"\t"$4"\t"$5}' >ENSG00000186092.bedsamtools depth -b ENSG00000186092.bed SRR1039510.Hisat_aln.sorted.bam >ENSG00000186092.bed.depth # 如何找到多比对的reads，flag值的理解# (0x100) 代表着多比对情况，所以直接用samtools view -f 0x100可以提取 multiple比对的 情况