我用的是virtualbox+ubuntu
记得修改内存、显存、存储、网络、共享文件夹
安装好系统后记得安装增强功能，开无缝模式
找到加载的磁盘，运行这个文件VBoxLinuxAdditions.run
下载SRA文件，sratools，现在没有Aspera这个东西了。
不过下载和解压都很快了。
prefetch SRR -p
fasterq-dump
然后用fastp质控
fastp -i SRR14213156_1.fastq -I SRR14213156_2.fastq -o out156_r1.fastq.gz -O out156_r2.fastq.gz -h h156.html -j j156.json -c -q 20 -n 7 -w 14 -z 6
然后就是比对，去hisat2官网找组装好的基因组，我用的这个包含转录组注释
hisat2 -p 6 --summary-file s167.txt -x ./grch38_tran/genome_tran -1 out167_r1.fastq.gz -2 out167_r2.fastq.gz -S align167.sam
samtools view -1 -@ 6 align165.sam -o align165.bam
samtools sort -@ 6 align165.bam -o sorted165.bam
fastq是序列数据
sam是比对数据，体积很大，转成二进制bam
然后是组装（-e模式）
for file in ./*.bam
do
stringtie -o $file.gtf -e -G GRCh38.gtf -p 6 $file
done
不同来源的gtf是不一样的，要选择对应版本，染色体名称不一样等等？
然后是提取count
ls *.bam.gtf | sed -e 's/sorted//; s/.bam.gtf//' | awk '{printf "sample%s\tsorted%s.bam.gtf\n",$1,$1}' > path1.txt
#sample17 sorted17.bam.gtf
#sample18 sorted18.bam.gtf
#sample19 sorted19.bam.gtf
#sample20 sorted20.bam.gtf
#sample21 sorted21.bam.gtf
#sample22 sorted22.bam.gtf
#sample23 sorted23.bam.gtf
提取tpm和fpkm
for file in ./*.bam.gtf
do
awk -F'"' '$1~/transcript/ {print $0}' $file | awk -F'"' '$5~/ref_gene/ {printf "%s\t%s\t%s\t%s\n",$4,$6,$10,$12}' > $file.txt
done
参考基因文件下载：
保姆级参考基因组及其注释下载方法（图文详解） - 知乎 (zhihu.com)
常用参考基因组—下载站点 - 简书 (jianshu.com)
关于参考基因组和注释 - 简书 (jianshu.com)
尝试了一下lncRNA鉴定，https://mp.weixin.qq.com/s/-q-415v4dm7BN_HLZ_jc7A
很重要的一点，在预测lncRNA之前用stringtie的--merge模式把序列合称为一个文件用于分析。
cd ~/Documents
#assembl first assembl against transcriptome
stringtie -o assembl.gtf -G GRCh38.gtf -p 6 sorted$1.bam
#compare with transcriptome
gffcompare -R -r GRCh38.gtf -o compare.gtf assembl.gtf
#summary for assemble class
awk '$3!~/class/ {print $3}' compare.gtf.assembl.gtf.tmap | sort -V | uniq -c
#class filter
awk '{if ($3=="u" || $3=="x" || $3=="i" || $3=="j" || $3=="o"){print $0}}' compare.gtf.assembl.gtf.tmap > filter1_byuxijo.tmap
#len and number of exon filter
awk '($6>1 && $10>=200){print $0}' filter1_byuxijo.tmap > filter_bylen.tmap
#get id
awk '{print $5}' filter_bylen.tmap > ID
#
grep -w -Ff ID -w assembl.gtf > filter.gtf
#exon gtf
awk '($3 == "exon") {print $0}' filter.gtf > filter_exon.gtf
#get exon fa based on len and uxijo filter against genome fasta
gffread -w exon.fa -g primary_assembly.fa filter_exon.gtf
#lncfinder pridiction
Rscript lf.R exon.fa
#cpat strange stdin and stdout
#cpat.py -x Human_Hexamer.tsv -d Human_logitModel.RData -g exon.fa -o out_cpat > cpat.log 2>&1
#ccp2 not with .txt
./CPC2/bin/CPC2.py -i exon.fa -o out_cpc2
#plek
python ./PLEK.1.2/PLEK.py -fasta exon.fa -out out_plek.txt -thread 10
#get ID_nc
grep -w 'NonCoding' out_lncf.txt > lncf_id.csv
awk -F '"' '{print $2}' lncf_id.csv > ID_lncf
#cp out_cpat.no_ORF.txt ID_cpat
awk '$8~/noncoding/ {print $1}' out_cpc2.txt > ID_cpc2
sed -n '/Non/s/>//p' out_plek.txt | awk '{print $3}' > ID_plek
#merge all ID files
cat ID_* | sort | uniq -c | awk '{if($1==3){print $2}}' > $1ID.last
#get exon.fa again based on non-coding
#grep -w -Ff ID.last -w assembl.gtf > filter2.gtf
#awk '($3 == "exon") {print $0}' filter2.gtf > filter_exon2.gtf
#gffread -w exon2.fa -g primary_assembly.fa filter_exon2.gtf
#transeq exon2.fa $1_protein.fa -frame 6
#./PfamScan/pfam_scan.pl -fasta protein.fa -dir ./PfamScan/ -outfile filter_bypfam$1.out
#grep -c '^>' exon.fa
#Rscript lncfinder.R exon.fa
#grep -w 'NonCoding' result.csv > lncfinder.csv
主要参考
http://journals.im.ac.cn/html/cjbcn/2015/9/gc15091271.htm
https://cloud.tencent.com/developer/article/1842208
网友评论